Skip to content
Snippets Groups Projects
Commit 1246c5be authored by weihe's avatar weihe
Browse files

revised errors on time mismatch between fluxes and footprints

parent ecf0a752
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
da/stilt/stilt.co2.simu.2015.Sep.bc.r
+ 121
− 88
View file @ 1246c5be
#------------------------------------------------------------------------------------------------
#i------------------------------------------------------------------------------------------------
# CO2 concentration simulation based on WRF-STILT/HYSPLIT-NAM12 footprints and
# SiB3/4 biosphere fluxes, CarbonTracker background fluxes and boundary conditions
#
......@@ -70,7 +70,7 @@ library(futile.logger)
logfile = paste("stilt_",procs,".log", sep="")
flog.appender(appender.file(logfile), name='logger.b')
flog.info("Launch WRF-STILT-based forward simulations %s", "file", name='logger.b')
flog.info("======Launch WRF-STILT-based forward simulations======", name='logger.b')
flog.info("The biosphere flux is from %s", bioflux_flag, name='logger.b')
flog.info("The boundary is from %s", boundary_flag, name='logger.b')
......@@ -204,7 +204,7 @@ newfile=T #output results into a new file
if(boundary_flag == "EMP")
{
bouf=paste(bounddir,y1,"/",m1,"/","CO2.v201209_v1.init.stilt.",y1,m1,".txt",sep="")
print(bouf)
#print(bouf)
val<-as.matrix(read.table(bouf,header=TRUE))
}
......@@ -221,6 +221,8 @@ for(mm in 1:length(pfbfns))
footnc=nc_open(fn)
foot=ncvar_get(footnc,"foot1",start=c(41,12,1),count=c(ncol2,nrow2,-1))
#print(mean(foot))
# get particles
#part=footp$partfoot
#partna=footp$partfootnames
......@@ -255,6 +257,8 @@ for(mm in 1:length(pfbfns))
# get the time stamp for each foot step, going back time given by "endtime"
xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
print("XXXxx?")
print(xx)
yy=xx+(-foottimes*3600)
cyy=as.character(yy)
yrlist=substring(cyy,1,4)
......@@ -268,7 +272,9 @@ for(mm in 1:length(pfbfns))
udaystring=unique(daystring)
yrmonstring=paste(yrlist,monlist, sep="")
uyrmonstring=unique(yrmonstring)
print("XXXudaystring 1st from ident name")
print(udaystring)
#current month
sibyr=substring(uyrmonstring[1],1,4)
sibmon=substring(uyrmonstring[1],5,6)
......@@ -292,24 +298,26 @@ for(mm in 1:length(pfbfns))
hgtarr=endpts[,5] + endpts[,6] #agl+grdht
# analyze the dates to recognize how many data to read
dd=as.character(endptsdate)
yrlist=substring(dd,1,4)
monlist=substring(dd,6,7)
daylist=substring(dd,9,10)
hrlist=substring(dd,12,13)
milist=substring(dd,15,16)
#dd=as.character(endptsdate)
#yrlist=substring(dd,1,4)
#monlist=substring(dd,6,7)
#daylist=substring(dd,9,10)
#hrlist=substring(dd,12,13)
# get unique months and days
daystring=paste(yrlist,monlist,daylist, sep="")
udaystring=unique(daystring)
yrmonstring=paste(yrlist,monlist, sep="")
uyrmonstring=unique(yrmonstring)
#daystring=paste(yrlist,monlist,daylist, sep="")
#udaystring=unique(daystring)
#yrmonstring=paste(yrlist,monlist, sep="")
#uyrmonstring=unique(yrmonstring)
#print("daylist=") #it's wrong to use daylist independently, but to use daystring (yyyymmdd)
#print(daylist)
#--------------------------------------------------------------
# 2-1: Read fluxes and boundary data
#--------------------------------------------------------------
ndays=length(udaystring)
if(boundary_flag=="CTNA")
bouarr=array(0,dim=c(ndays,120,90,25,8)) #boundary : lon,lat,hgt,time
......@@ -454,9 +462,9 @@ for(mm in 1:length(pfbfns))
inxd=ceiling(hh/24)
inxh=ceiling( (hh-(inxd-1)*24)/3 ) #every three hours the same
print(hh)
print(inxh)
print(dim(ocn))
#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]
......@@ -475,28 +483,7 @@ for(mm in 1:length(pfbfns))
gppflux[,,]=replace(gppflux[,,],is.na(gppflux[,,]),0)
recflux[,,]=replace(recflux[,,],is.na(recflux[,,]),0)
#for(hh in ntimes2:1)
#{
# dateall2[ntimes2-hh+1]=time2[hh]
# inxd=ceiling(hh/8)
# inxh=hh-(inxd-1)*8
# ocnflux[,,ntimes2-hh+1]=ocn[inxd,52:117,113:152,inxh] #transform matrix, delete this
# fosflux[,,ntimes2-hh+1]=fos[inxd,52:117,113:152,inxh]
# firflux[,,ntimes2-hh+1]=fir[inxd,52:117,113:152,inxh]
# if(bioflux_flag == "SiBCASA")
# {
# gppflux[,,ntimes2-hh+1]=gpp[inxd,52:117,113:152,inxh]
# recflux[,,ntimes2-hh+1]=rec[inxd,52:117,113:152,inxh]
# }
# if(bioflux_flag == "CT_OPT")
# bioflux[,,ntimes2-hh+1]=bio[inxd,52:117,113:152,inxh]
#}
# replace NA values with 0, as NA was used as zero values for datasets
ocnflux[,,]=replace(ocnflux[,,],is.na(ocnflux[,,]),0)
fosflux[,,]=replace(fosflux[,,],is.na(fosflux[,,]),0)
firflux[,,]=replace(firflux[,,],is.na(firflux[,,]),0)
......@@ -514,8 +501,6 @@ for(mm in 1:length(pfbfns))
if(bioflux_flag == "CT_OPT")
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')
#--------------------------------------------------------------
......@@ -525,8 +510,7 @@ for(mm in 1:length(pfbfns))
#**********************************
# biosphere fluxes convolving
#**********************************
# hourly flux
datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600
# hourly fxx-0.5*24*3600-ateall[1]+0.5*3600-(0:(ntimes-1))*3600
ixflux=(1:ntimes)[(datevalid<=xx)][1:endtime] #time backward
if(bioflux_flag == "SiB3")
......@@ -549,8 +533,11 @@ for(mm in 1:length(pfbfns))
# involve scaling factors to neeflux[,,ixflux], partly, for example as follow
# xx = "2010-01-24 18:10:00 UTC" || b1 = "2010-01-19 12:00:00 UTC", b3 = "2010-01-29 12:00:00 UTC"
xxleft=xx-ndays*24*3600
#xxleft=xx-ndays*24*3600
xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
yy=xx+(-foottimes*3600)
xxleft=min(yy)
#flog.info("xxleft = %s, b2=%s, b2-xxleft=%s", xxleft, b2, b2-xxleft,name='logger.b')
......@@ -651,6 +638,11 @@ for(mm in 1:length(pfbfns))
fosout=sum(fostmp,na.rm=TRUE)
firout=sum(firtmp,na.rm=TRUE)
print("XXX Line 654")
print(mean(foot))
print(mean(fosflux))
print(fosout)
if(bioflux_flag == "CT_OPT")
{
biotmp=foot*bioflux[,,ixflux]*1e6
......@@ -673,60 +665,101 @@ for(mm in 1:length(pfbfns))
lonmax = 180 #-63.5
npts=dim(endpts)[1]
# calulate concentrations
# calculate concentrations
pbou=0
nval=0
dmin=min(as.integer(substr(udaystring,7,8)))
#dmin=min(as.integer(substr(aystring,7,8)))
#dmax=max(as.integer(substr(udaystring,7,8)))
dmin=min(endptsdate)
dmax=max(endptsdate)
xx=ISOdatetime(inityr,initmo,initdy,inithr,initmi,0,tz="UTC")
yy=xx+(-foottimes*3600)
print("XXXudaystring for selected flux period")
print(udaystring)
print("XXXdmin particles")
print(dmin)
print("XXXdmax particles")
print(dmax)
print("XXXdobs time")
print(xx)
dd1=as.character(endptsdate)
yrlist=substring(dd1,1,4)
monlist=substring(dd1,6,7)
daylist=substring(dd1,9,10)
hrlist=substring(dd1,12,13)
milist=substring(dd1,15,16)
for(p in 1:npts)
{
dy=0
goal=as.integer(daylist[p])
for(uu in 1:ndays)
#dy=0
#goal=as.integer(daylist[p])
dy=ceiling( (as.integer(difftime(xx,endptsdate[p],units='hours')) )/24 )+1
#for(uu in 1:ndays)
#{
# dstr=udaystring[uu]
# ddy=as.integer(substr(dstr,7,8))
#print(ddy)
#print(goal)
# if(ddy<=goal)
# dy=dy+1
#}
#print("XXXdy=")
#print(dy)
if(boundary_flag == "CTNA")
{
dstr=udaystring[uu]
ddy=substr(dstr,7,8)
if(ddy<=goal)
dy=dy+1
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)
}
if(boundary_flag == "CTNA")
{
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
for(l in 1:25)
{
if(k0 > levelhgt[l] && k0 <= levelhgt[25])
k=l+1
if(k0 > levelhgt[25])
k=25
}
# height matching
k0=hgtarr[p]
k=1
for(l in 1:25)
{
if(k0 > levelhgt[l] && k0 <= levelhgt[25])
k=l+1
if(k0 > levelhgt[25])
k=25
}
# 3-hourly matching /agt/alt
hpos=ceiling((as.integer(hrlist[p])+0.000001+as.integer(milist[p])/60.0 )/3.0)
# 3-hourly matching /agt/alt
hpos=ceiling((as.integer(hrlist[p])+1e-6+as.integer(milist[p])/60.0 )/3.0)
tt=bouarr[dy,i,j,k,hpos]
print("dy")
print(dy)
print("hpos")
print(hpos)
tt=bouarr[dy,i,j,k,hpos]
if(length(tt)==1) #why sometimes we get a unnormal array?
{
pbou=pbou+tt #sum
nval=nval+1
}
if(length(tt)==1) #why sometimes we get a unnormal array?
{
pbou=pbou+tt #sum
nval=nval+1
}
}
fbou=pbou/nval
}
print("pbou")
print(pbou)
print(nval)
fbou=pbou/nval
}
print("XXXfbou")
print(fbou)
#case 2: data from Arlyn's emperical curtain
if(boundary_flag == "EMP")
{
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment