Commit a91db84d authored by weihe's avatar weihe
Browse files

to involve 2011 January or more footprints

parent e4091cc4
......@@ -25,7 +25,7 @@ procs= as.integer(args[1])
#inputs and outputs
stiltfile = paste("stilt_",procs,".rc", sep="")
curdir="/Storage/CO2/wei/ctdas-stilt-parallel-f-only-newcov-bugfix-updatedfoot/exec/da/stilt/"
curdir="/Storage/CO2/wei/ctdas-stilt-parallel-f-only-newcov-bugfix-updatedfoot-sumwinmdm-bcCTE/exec/da/stilt/"
fn=paste(curdir, stiltfile,sep="")
conf=as.matrix(read.table(fn,header=FALSE))
......@@ -37,8 +37,8 @@ samdir = conf[5,3]
outdir = conf[6,3]
# data choice flag
bioflux_flag = "SiB3" #"SiB3", "SiBCASA", "CT_OPT"
boundary_flag= "CT" #"CT", "EMP"
bioflux_flag = "SiBCASA" #"SiB3", "SiBCASA", "CT_OPT"
boundary_flag= "CTE" #"CTNA", "CTE", "EMP"
# optimization option (flux only, or flux + boundary)
opt_para_flag = "FLUX" # "FLUX_BOU","FLUX"
......@@ -133,19 +133,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")
......@@ -182,13 +188,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)
......@@ -282,7 +300,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
......@@ -323,9 +345,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")
......@@ -352,13 +377,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"))
}
#--------------------------------------------------------------
......@@ -479,6 +513,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')
#--------------------------------------------------------------
......@@ -628,33 +663,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
......@@ -677,8 +723,23 @@ for(mm in 1:length(pfbfns))
nval=nval+1
}
}
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
}
}
fbou=pbou/nval
#################################################
# 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