Commit 906c2960 authored by weihe's avatar weihe
Browse files

time-involved BC sensitivity for WRF-STILT

parent 7186f313
#------------------------------------------------------------------------------------------------
# CO2 concentration simulation based on WRF-STILT/HYSPLIT-NAM12 footprints and
# SiB3/4 biosphere fluxes, CarbonTracker background fluxes and boundary conditions
#
# created by W.He on Feb 20, 2015
#------------------------------------------------------------------------------------------------
source("/Users/wei/co2simu/Rcode/lpdm/rsource/load.ncdf4.r")
source("/Users/wei/co2simu/Rcode/empirical_boundary/src/id2info.r")
source("/Users/wei/co2simu/Rcode/empirical_boundary/src/julian.r")
source("/Users/wei/co2simu/Rcode/rsource/time.r")
source("/Users/wei/co2simu/Rcode/rsource/ddate.r")
source("/Users/wei/co2simu/Rcode/rsource/datelist.r")
source("/Users/wei/co2simu/Rcode/stilt-cvs-20090417/stiltR/sourceall.r")
source("/Users/wei/co2simu/Rcode/stilt-cvs-20090417/stiltR/Trajecfoot.r")
source("/Users/wei/co2simu/Rcode/stilt-cvs-20090417/stiltR/getgridp.r")
source("/Users/wei/co2simu/Rcode/rsource/assignr.r")
source("/Users/wei/co2simu/Rcode/output.ncdf4.v2.2.r")
args <- commandArgs(trailingOnly = TRUE)
procs= as.integer(args[1])
#inputs and outputs
stiltfile = paste("stilt_",procs,".rc", sep="")
curdir="/Storage/CO2/wei/ctdas-stilt-BCfp-randominit-alt3km-time/exec/da/stilt/"
fn=paste(curdir, stiltfile,sep="")
conf=as.matrix(read.table(fn,header=FALSE))
sibdir = conf[1,3]
bgfdir = conf[2,3]
footpdir= conf[3,3]
bounddir= conf[4,3]
samdir = conf[5,3]
outdir = conf[6,3]
# data choice flag
bioflux_flag = "SiBCASA" #"SiB3", "SiBCASA", "CT_OPT"
boundary_flag= "CTNA" #"CTNA", "CTE", "EMP"
# optimization option (flux only, or flux + boundary)
opt_para_flag = "FLUX" # "FLUX_BOU","FLUX"
#----------------------------------------------------------------------------------------------
#Convolve footprints with sib hourly fluxes, for observations from both Aircraft and towers
#----------------------------------------------------------------------------------------------
endtime=240 #hourly for 10 days
foottimes=seq(0,endtime,1) #vector of times (backtimes) in hours between which footprint is computed
zbot=0 #lower vertical bound for influence projection, in meters agl
ztop=0 #upper vertical bound for influence projection, in meters agl
#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
ResLon2=1 #resolution in degrees longitude
ResLat2=1 #resolution in degrees latitude
# constants for level-height transfrom
levelhgt=c(34.5,111.9,256.9,490.4,826.4,1274.1,1839.0,2524.0,3329.9,4255.6,5298.5,6453.8,7715.4,9076.6,
+ 10533.3,12108.3,13874.2,15860.1,18093.2,20590.0,24247.3,29859.6,35695.0,42551.5,80000.0)
#data from http://www.esrl.noaa.gov/gmd/ccgg/carbontracker-ch4/documentation_tm5.html
TBegin <- proc.time()
library(ncdf4)
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======", 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')
nsam=as.numeric(conf[7,3])
ndayscycle=as.numeric(conf[8,3])
startdate=conf[9,3]
year=as.numeric(substring(startdate,1,4))
month=as.numeric(substring(startdate,6,7))
day=as.numeric(substring(startdate,9,10))
b0 = ISOdatetime(year,month,day,0,0,0,tz="UTC")
startdate=conf[10,3]
year=as.numeric(substring(startdate,1,4))
month=as.numeric(substring(startdate,6,7))
day=as.numeric(substring(startdate,9,10))
if(ndayscycle!=10)
flog.fatal("Sorry, only support the case nndayscycle equals 10", name='logger.b')
b2 = ISOdatetime(year,month,day,0,0,0,tz="UTC")
b1 = b2 - ndayscycle*24*3600 # for footprint
b3 = b1 + 2*ndayscycle*24*3600 # for flux
tb1 = paste(substring(b1,1,4),substring(b1,6,7),substring(b1,9,10),sep="")
tb2 = paste(substring(b2,1,4),substring(b2,6,7),substring(b2,9,10),sep="")
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))
# make CTDAS to generate domain for North America, read data partly?
flog.info("Reading scaling factor files", name='logger.b')
for(i in 0:(nsam-1)) #parameters.000.2010010100_2010011100.nc
{
if (i<10)
ii = paste("00",i,sep="")
if (i<100 && i>=10)
ii = paste("0",i,sep="")
if (i>=100)
ii = i
if (b1<b0) #b1==b0, revise on Aug 5,2015
{
scalefacarr1[,,i+1] = 1
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,"parametervalues_bc")
scalefacarr_bc1[i+1] = scalefac
nc_close(ncf)
}
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,"parametervalues_bc")
scalefacarr_bc2[i+1] = scalefac
nc_close(ncf)
}
#---------------------------------------------------------------------------------
# Centering at the state vector especially the ready-optimized cycle, look for data
# 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(yy in 1:length(uyr))
{
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")
obs_id <- ncvar_get(ncf,"obs_id")
# filter before loop
pfbfns=NULL
fn=NULL
eventids=NULL
obsids=NULL
for(mm in 1:length(fns))
{
eid=as.numeric(substring(fns[mm],48,53))
for(i in 1:length(eventidarr))
{
if(eid==eventidarr[i])
{
pfbfns = c(pfbfns,fns[mm])
obsids = c(obsids,obs_id[i])
eventids = c(eventids,eventidarr[i])
break
}
}
}
nc_close(ncf)
#evid matching result
n_obs=length(eventidarr)
n_foot=length(fns)
n_used=length(pfbfns)
log=paste("Number of event ids from ctdas: ",n_obs,", Number of matched: ",n_used)
flog.info(log, name='logger.b') #June 23
log=paste("For state vectors from ",b1,"to",b2,",",length(pfbfns),"observations have been found")
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
#--------------------------------------------
# flog.info('obsids: %s',obsids[mm], name='logger.b')
# flog.info('eventids: %s',eventids[mm], name='logger.b') #June 23
yr=substring(pfbfns[mm],6,9)
mo=substring(pfbfns[mm],11,12)
fn=paste(footpdir,yr,"/",mo,"/",pfbfns[mm],sep="")
footp=load.ncdf(fn)
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
#colnames(part)=c("time","index","lat","lon","agl","grdht","foot","temp","temp0","swrad","zi","dens","pres","dmass") #ndmass
# get info form path srings
ident=substring(pfbfns[mm],6,46)
eventid=substring(pfbfns[mm],48,53)
print(ident)
flog.info("Convolving for %s",ident, name='logger.b')
#info=id2info(ident)
#print(info)
#foot1=Trajecfoot(ident=ident,part=part, pathname="",foottimes=foottimes,zlim=c(zbot,ztop),fluxweighting=NULL,coarse=1,vegpath=vegpath,numpix.x=ncol2,numpix.y=nrow2,lon.ll=LLLon2,lat.ll=LLLat2,lon.res=ResLon2,lat.res=ResLat2)
#foot1 40 66 240
#foot=array(NA, dim=c(66,40,240))
#for (i in 1:240)
# foot[,,i]=t(as.matrix(foot1[,,i])) #113:152
nc_close(footnc)
if(length(foot)>100)
{
inityr=as.numeric(substring(ident,1,4))
initmo=as.numeric(substring(ident,6,7))
initdy=as.numeric(substring(ident,9,10))
inithr=as.numeric(substring(ident,12,13))
initmi=as.numeric(substring(ident,15,16))
inihgt=as.numeric(substring(ident,39,41))
# 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) #reversed order
cyy=as.character(yy)
yrlist=substring(cyy,1,4)
monlist=substring(cyy,6,7)
daylist=substring(cyy,9,10)
hrlist=substring(cyy,12,13)
milist=substring(cyy,15,16)
# get unique months and days
daystring=paste(yrlist,monlist,daylist, sep="")
udaystring=unique(daystring)
yrmonstring=paste(yrlist,monlist, sep="")
uyrmonstring=unique(yrmonstring)
udaystring=rev(udaystring)
print("XXXudaystring 1st from ident name")
print(udaystring)
#current month
sibyr=substring(uyrmonstring[1],1,4)
sibmon=substring(uyrmonstring[1],5,6)
#----------------------------------------------------------------------------------
# STEP 2: Read boundary conditions & use end points tracing to
# 1) get responding fluxes & convolve with footprints
# 2) get the "actural" concentrations
#----------------------------------------------------------------------------------
# get endpoints
endpts=footp$endpts
endptna=footp$endptsnames
#endpts=ncvar_get(footnc,"endpts")
#endptna=ncvar_get(footnc,"endptsnames")
colnames(endpts)=endptna #=c("time","index","lat","lon","agl","grdht","temp","pres")
endptsdate=footp$endptsdate #2010-08-04 20:18:00 UTC
#endptsdate=ncvar_get(footnc,"endptsdate")
latarr=endpts[,3]
lonarr=endpts[,4]
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)
# get unique months and days
#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
if(boundary_flag=="CTE")
bouarr=array(0,dim=c(ndays,360,180,25,8))
#biospheric fluxes
if(bioflux_flag == "SiB3")
{
nrow=181 #181
ncol=360 #361 #288
gpp=array(NA, dim=c(ndays,ncol,nrow,24))
rec=array(NA, dim=c(ndays,ncol,nrow,24))
}
#other fluxes #(360,180,8)
nrow=180
ncol=360
ocn=array(NA, dim=c(ndays,ncol,nrow,8))
fos=array(NA, dim=c(ndays,ncol,nrow,8))
fir=array(NA, dim=c(ndays,ncol,nrow,8))
if(bioflux_flag == "SiBCASA")
{
gpp=array(NA, dim=c(ndays,ncol,nrow,8))
rec=array(NA, dim=c(ndays,ncol,nrow,8))
}
if(bioflux_flag == "CT_OPT")
bio=array(NA, dim=c(ndays,ncol,nrow,8))
ntimes=ndays*24
for(d in 1:ndays)
{
datestr=udaystring[d]
yr=substr(datestr,1,4)
mn=substr(datestr,5,6)
dy=substr(datestr,7,8)
print("#####################")
print(datestr)
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")
{
biof=load.ncdf(paste(sibdir,"SiB3.hourly.flux1x1.global.",yr,mn,dy,".nc",sep=""))
gpp[d,,,]=biof$gpp
rec[d,,,]=biof$rtotal
#pla[d,,,]=biof$ocs.gpp
#soi[d,,,]=biof$ocs.soil
remove(list=c("biof"))
}
if(bioflux_flag == "SiBCASA")
{
biof=load.ncdf(paste(sibdir,"SiBCASA.3hourly.flux1x1.global.",yr,mn,dy,".nc",sep=""))
gpp[d,,,]=biof$gpp
rec[d,,,]=biof$resp
remove(list=c("biof"))
#flog.info('d: %s, yr: %s, mn: %s, dy: %s',d,yr,mn,dy,name='logger.b')
}
#flog.info("mean gpp= %f", mean(gpp[d,,,]), name='logger.b')
if(bioflux_flag == "CT_OPT")
bio[d,,,]=bgf$bio.flux.opt
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("bgf"))
}
#--------------------------------------------------------------
# 2-2: prepare data for calculations
#--------------------------------------------------------------
dateall=rep(ISOdatetime(0,0,0,0,0,0,tz="UTC"), ntimes)
if(bioflux_flag == "SiB3" || bioflux_flag == "SiBCASA")
{
gppflux=array(NA, dim=c(ncol2,nrow2,ntimes))
recflux=array(NA, dim=c(ncol2,nrow2,ntimes))
}
if(bioflux_flag == "CT_OPT")
bioflux=array(NA, dim=c(ncol2,nrow2,ntimes))
ocnflux = array(NA, dim=c(ncol2,nrow2,ntimes))
fosflux = array(NA, dim=c(ncol2,nrow2,ntimes))
firflux = array(NA, dim=c(ncol2,nrow2,ntimes))
neeflux = array(NA, dim=c(ncol2,nrow2,ntimes))
neeflux1 = array(NA, dim=c(ncol2,nrow2,ntimes))
neefluxarr= array(NA, dim=c(ncol2,nrow2,ntimes,nsam))
bc1 = array(NA, dim=c(ntimes))
bcarr= array(NA, dim=c(ntimes,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))
#----------------------------------------------------
yr=rep(sibyr,ntimes)
mon=rep(sibmon,ntimes)
hrs=seq(1,ntimes,1)
dy=ceiling(hrs/24)
hr=(hrs-(dy-1)*24)*1-1
time=ISOdatetime(yr,mon,dy,hr,0,0,tz="UTC")
for(hh in ntimes:1)
{
dateall[ntimes-hh+1]=time[hh]
if(bioflux_flag == "SiB3")
{
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]
}
inxd=ceiling(hh/24)
inxh=ceiling( (hh-(inxd-1)*24)/3 ) #every three hours the same
#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]
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]
#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]
}
# replace NA values with 0, as NA was used as zero values for datasets
gppflux[,,]=replace(gppflux[,,],is.na(gppflux[,,]),0)
recflux[,,]=replace(recflux[,,],is.na(recflux[,,]),0)
ocnflux[,,]=replace(ocnflux[,,],is.na(ocnflux[,,]),0)
fosflux[,,]=replace(fosflux[,,],is.na(fosflux[,,]),0)
firflux[,,]=replace(firflux[,,],is.na(firflux[,,]),0)
if(bioflux_flag == "CT_OPT")
bioflux[,,]=replace(bioflux[,,],is.na(bioflux[,,]),0)
if(bioflux_flag == "SiB3")
neeflux[,,]=recflux[,,]-gppflux[,,] # nee*1e6 for SiBCASA?
if(bioflux_flag == "SiBCASA")
neeflux[,,]=(recflux[,,]-gppflux[,,])*1e6
if(bioflux_flag == "CT_OPT")
neeflux[,,]= bioflux[,,]*1e6 # for scaling CT flux
#flog.info("nee flux, max=%f, min=%f, mean=%f", max(neeflux),min(neeflux),mean(neeflux), name='logger.b')
#--------------------------------------------------------------
# 2-3: Convolving footprints with fluxes to get delta co2
#--------------------------------------------------------------
#**********************************
# biosphere fluxes convolving
#**********************************
# hourly flux
#print("XXXdateall")
#print(dateall) #reversed order
# # note: ixflux should not be from 1:ntimes, should be more accurate in hours
# dd=as.character(min(yy))
# #yrlist=substring(dd,1,4)
# #monlist=substring(dd,6,7)
# #daylist=substring(dd,9,10)
# hr=substring(dd,12,13)
# mi=substring(dd,15,16)
# n = floor(as.numeric(hr) + as.numeric(mi)/60)
# n = 24-n
# print("XXXn")
# print(n)
# datevalid=dateall[1]+0.5*3600-(0:(ntimes-1))*3600
# #ixflux=(1:ntimes)[(datevalid<=xx)][1:endtime] #time backward
# #ixflux=(ntimes:1)[(datevalid<=xx)][(endtime+n-1):n]
# #ixflux=(ntimes:1)[(datevalid<=xx)][(endtime+n-1):(n)]
# ixflux<-c(n:endtime+n-1)
## ixflux=(1:ntimes)[(datevalid<=xx)][n:(endtime+n-1)]
#dd=as.character(min(yy))
dd=as.character(xx)
#yrlist=substring(dd,1,4)
#monlist=substring(dd,6,7)
#daylist=substring(dd,9,10)
hr=substring(dd,12,13)
mi=substring(dd,15,16)
#flog.info('dd: %s, length dd: %s,hr: %s, mi: %s',dd,length(dd),hr,mi,name='logger.b')
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
#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')
ixflux <-array(n:(endtime+n-1),dim=c(240))
#flog.info('ixflux endtime: %s,%s,%s',endtime,n,as.numeric(hr)+(as.numeric(hr)/60),name='logger.b')
#for (ttt in 1:240)
# flog.info('ixflux: %s,%s',ttt,ixflux[ttt],name='logger.b')
if(bioflux_flag == "SiB3") #footprint is backward
{
gpptmp=foot*gppflux