Skip to content
Snippets Groups Projects
Commit 9beec836 authored by Peters, Wouter's avatar Peters, Wouter
Browse files

The templates for jobs and rc-files used by start_ctdas.sh are now in a... 

The templates for jobs and rc-files used by start_ctdas.sh are now in a separate templates folder. This allows each user to modify their templates for running if needed, and to add templates for their own purpose.
parent d498349c
Branches
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
stilt.nc4.sib.co2.Feb.2015.May13.r deleted 100755 → 0
+ 0
− 600
View file @ d498349c
# CarbonTracker Data Assimilation Shell (CTDAS) Copyright (C) 2017 Wouter Peters.
# Users are recommended to contact the developers (wouter.peters@wur.nl) to receive
# updates of the code. See also: http://www.carbontracker.eu.
#
# This program is free software: you can redistribute it and/or modify it under the
# terms of the GNU General Public License as published by the Free Software Foundation,
# version 3. This program is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
# FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License along with this
# program. If not, see <http://www.gnu.org/licenses/>.
#------------------------------------------------------------------------------------------------
# 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.r")
# #inputs and outputs
# sibdir="/Users/weihe/Documents/Concensimulation/sib3_fluxes/day1×1/"
# bgfdir="/Users/weihe/Documents/Concensimulation/bgfluxes/3-hourly/"
# footpdir="/Users/weihe/Documents/Concensimulation/stilt_footprints/" #read arlyn's calculated data
# bounddir="/Users/weihe/Documents/Concensimulation/co2_boundary/co2_total/"
# outdir="/Users/weihe/Documents/Concensimulation/resultsforctdas/"
#
# #scaling factors for 10-days time scale
# samdir="/Users/weihe/Documents/example_ctdas/input/"
stiltfile="stilt.rc"
curdir="/Storage/CO2/wei/ctdas-stilt-proto/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]
outdir = conf[5,3] #output
samdir = conf[6,3] #input
#----------------------------------------------------------------------------------------------
#Convolve footprints with sib hourly fluxes, for observations from both Aircraft and towers
#----------------------------------------------------------------------------------------------
endtime=240
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)
#flog.appender(appender.file("/Users/weihe/Documents/Concensimulation/resultstest/stilt.log"), name='logger.b')
#flog.info("This writes to a logging %s", "file", name='logger.b')
#flog.warn("This statement has higher severity", name='logger.b')
#flog.fatal("This one is really scary", name='logger.b')
#--------------------------------------------------------------------------------------------------------
# site level: read scaling factors for all sites once
# member level: read boundary for all ensemble members once
#nsam= 150
#ndayscycle=10 # one cycle of footprints
#year=2010
#month=1
#day=19
nsam=as.numeric(conf[7,3]) #get from stilt.txt file =>replace with a dictionary style, not id order; also check parameters and report status
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))
#if(ndayscycle!=10)
# flog.warn("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))
# 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
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)) #temporary 42 #real52:117,113:152,start=c(52,113),count=c(66,40)
scalefacarr1[,,i+1]=scalefac
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
}
#---------------------------------------------------------------------------------
# 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
m1=substring(b2,6,7)
m2=substring(b3,6,7)
umon=unique(c(m1,m2))
#flog.info("Determing ready-use footprint files", name='logger.b')
fns=NULL
for(mm in 1:length(umon))
{
pfbpath=paste(footpdir,year,"/",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.***.nc files(e.g, sample_coordinates_2010011100_2010012100.nc)
ncf <- nc_open(paste(samdir,"sample_coordinates_",tb2,"00","_",tb3,"00",".nc",sep=""))
eventidarr <- ncvar_get(ncf,"obs_num")
# filter before loop
pfbfns=NULL
fn=NULL
for(mm in 1:length(fns))
{
#select: Uniform time to make it start from the same point
#yr=substring(fns[mm],6,9) #"2010-01-09 09:00:00 UTC"
#mo=substring(fns[mm],11,12)
#dy=substring(fns[mm],14,15)
#bISO = ISOdate(yr,mo,dy)
#diff1=as.numeric(bISO - b2)
#diff2=as.numeric(bISO - b3)
#if (diff1 >= 0 && diff2 <= 0) # in the time span
# pfbfns = c(pfbfns,fns[mm])
eid=as.numeric(substring(fns[mm],48,53))
for(i in 1:length(eventidarr))
{
if(eid==eventidarr[i])
{
pfbfns = c(pfbfns,fns[mm])
break
}
}
}
#log=paste("For state vectors from ",b1,"to",b2,",",length(pfbfns),"observations have been found")
#flog.info(log, name='logger.b')
# for fluxes, the current 10 days using scaling factors
#flog.info("Start convolving for all footprint files", name='logger.b')
newfile=T #output results into a new file
for(mm in 1:length(pfbfns))
{
# Read pfb file
#mm=1
mo=substring(pfbfns[mm],11,12)
fn=paste(footpdir,year,"/",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))
# Read footprint directly
#ft=footp$foot1 #120*70*240
#foot=array(NA, dim=c(40,66,240))
#for (i in 1:240)
# foot[,,i]=t(as.matrix(ft[41:106,12:51,i])) #113:152
# 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) #interpret info from file names
#--------------------------------------------------------
# STEP 1: Creates footprint for individual particle runs
#--------------------------------------------------------
#foot1=Trajecfoot(ident=ident,part=part, pathname="",foottimes=foottimes,zlim=c(zbot,ztop),fluxweighting=NULL,coarse=1,vegpath=vegpath,
# numpix.x=ncol1,numpix.y=nrow1,lon.ll=LLLon1,lat.ll=LLLat1,lon.res=ResLon1,lat.res=ResLat1)
# for different spatial resolutions between biosphere fluxes and background fluxes, two footprints are planned to be used, and to the end, we picked up the concentration values.
#foot2=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)
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)) #by He
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")
yy=xx+(-foottimes*3600)
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)
#kk=length(uyrmonstring) #needs July (former month)
#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)
milist=substring(dd,15,16)
# get unique months and days
daystring=paste(yrlist,monlist,daylist, sep="")
udaystring=unique(daystring)
yrmonstring=paste(yrlist,monlist, sep="")
uyrmonstring=unique(yrmonstring)
#--------------------------------------------------------------
# 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
#biospheric fluxes
nrow_flux1=181
ncol_flux1=361 #288
gpp=array(NA, dim=c(ndays,ncol_flux1,nrow_flux1,24))
rec=array(NA, dim=c(ndays,ncol_flux1,nrow_flux1,24))
pla=array(NA, dim=c(ndays,ncol_flux1,nrow_flux1,24))
soi=array(NA, dim=c(ndays,ncol_flux1,nrow_flux1,24))
#other fluxes #(360,180,8)
nrow_flux2=180
ncol_flux2=360
ocn=array(NA, dim=c(ndays,ncol_flux2,nrow_flux2,8))
fos=array(NA, dim=c(ndays,ncol_flux2,nrow_flux2,8))
fir=array(NA, dim=c(ndays,ncol_flux2,nrow_flux2,8))
ntimes1=ndays*24
ntimes2=ndays*8
for(d in 1:ndays)
{
datestr=udaystring[d]
yr=substr(datestr,1,4)
mn=substr(datestr,5,6)
dy=substr(datestr,7,8)
#bou=load.ncdf(paste(bounddir,"CT2013B.molefrac_nam1x1_",yr,"-",mn,"-",dy,".nc",sep=""))
bou=load.ncdf(paste(bounddir,"CT2013B.molefrac_glb3x2_",yr,"-",mn,"-",dy,".nc",sep=""))
# co2(date, level, lat, lon) / ocn_flux_opt(date, lat, lon)
biof=load.ncdf(paste(sibdir,"SiB3.hourly.flux1x1.global.",yr,mn,dy,".nc",sep=""))
bgf=load.ncdf(paste(bgfdir,"CT2013B.flux1x1.",yr,mn,dy,".nc",sep=""))
bouarr[d,,,,]=bou$co2 # boundary CO2
gpp[d,,,]=biof$gpp
rec[d,,,]=biof$rtotal
pla[d,,,]=biof$ocs.gpp
soi[d,,,]=biof$ocs.soil
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","biof","bgf"))
}
#--------------------------------------------------------------
# 2-2: prepare data for calculations
#--------------------------------------------------------------
dateall1=rep(ISOdatetime(0,0,0,0,0,0,tz="UTC"), ntimes1)
dateall2=rep(ISOdatetime(0,0,0,0,0,0,tz="UTC"), ntimes2)
gppflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
recflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
plaflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
soiflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
neeflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
neeflux1=array(NA, dim=c(ncol2,nrow2,ntimes1))
#----------------------------------------------------
neefluxarr=array(NA, dim=c(ncol2,nrow2,ntimes1,nsam))
neeoutarr=array(NA, dim=c(nsam))
fsimuarr=array(NA, dim=c(nsam))
#----------------------------------------------------
ocnflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
fosflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
firflux=array(NA, dim=c(ncol2,nrow2,ntimes1))
yr=rep(sibyr,ntimes1)
mon=rep(sibmon,ntimes1)
hrs=seq(1,ntimes1,1)
dy=ceiling(hrs/24)
hr=(hrs-(dy-1)*24)*1-1
time1=ISOdatetime(yr,mon,dy,hr,0,0,tz="UTC")
yr=rep(sibyr,ntimes2)
mon=rep(sibmon,ntimes2)
hrs=seq(1,ntimes2,1)
dy=ceiling(hrs/8)
hr=(hrs-(dy-1)*8)*3-1
time2=ISOdatetime(yr,mon,dy,hr,0,0,tz="UTC")
for(hh in ntimes1:1)
{
dateall1[ntimes1-hh+1]=time1[hh]
inxd=ceiling(hh/24)
inxh=hh-(inxd-1)*24
gppflux[,,ntimes1-hh+1]=gpp[inxd,52:117,113:152,inxh] #transform matrix, 41:93; delete this on May 4,2015
recflux[,,ntimes1-hh+1]=rec[inxd,52:117,113:152,inxh]
plaflux[,,ntimes1-hh+1]=pla[inxd,52:117,113:152,inxh]
soiflux[,,ntimes1-hh+1]=soi[inxd,52:117,113:152,inxh]
}
neeflux[,,]=recflux[,,]-gppflux[,,]
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]
}
#--------------------------------------------------------------
# 2-3: Convolving footprints with fluxes to get delta co2
#--------------------------------------------------------------
# sib3 flux, 1-hourly
datevalid=dateall1[1]+0.5*3600-(0:(ntimes1-1))*3600 #need to change, ntimes=1488
ixflux=(1:ntimes1)[(datevalid<=xx)][1:endtime] #time backward
#dateall[ixflux]
#datevalid from dateall1, from time1, from SiB flux; xx from from footprint
#need large amount of memories
ocstmp=foot*plaflux[,,ixflux]
soiltmp=foot*soiflux[,,ixflux]
gpptmp=foot*gppflux[,,ixflux]
recotmp=foot*recflux[,,ixflux] #dim(neetmp): 40 53 240
#################################################
# determine which fluxes need scaling factors
#################################################
# 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
for (i in 1:nsam)
{
if(xxleft<b2) #flux with scaling factors
{
diff=as.numeric(difftime(b2,xxleft,units='hours')) # needs scaling for first half
diff=ceiling(diff)
for (hh in 1:(diff))
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr1[,,i]) #delte "t" transform
for (hh in (diff+1):ntimes1)
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
neefluxarr[,,,i] = neeflux1[,,]
}
else
{
for (hh in 1:ntimes1)
neeflux1[,,hh] = neeflux[,,hh]*(scalefacarr2[,,i])
neefluxarr[,,,i] = neeflux1[,,]
}
}
# REMEMBER: the spatial resolutios of the scaling factors and fluxes are different, whiched needs to be uniformed afterwards
# delta co2 on NEE for all ensemble members
for(i in 1:nsam)
{
neetmp=foot*neefluxarr[,,ixflux,i]
neeout=sum(neetmp,na.rm=TRUE)
neeoutarr[i]=neeout
}
#------------------------------------------------
ocsout=sum(ocstmp,na.rm=TRUE)
soilout=sum(soiltmp,na.rm=TRUE)
gppout=sum(gpptmp,na.rm=TRUE)
recoout=sum(recotmp,na.rm=TRUE)
#remove(list=c("gppflux","recflux","soiflux","plaflux","gpptmp","recotmp","soiltmp","ocstmp","dateall1"))
gc()
# bg fluxes, 3-hourly
datevalid=dateall2[1]+1.5*3600-(0:(ntimes2-1))*3*3600
ixflux=(1:ntimes2)[(datevalid<=xx)][1:endtime]
#need large amount of memories
ocntmp=foot*ocnflux[,,ixflux]*1e6
fostmp=foot*fosflux[,,ixflux]*1e6
firtmp=foot*firflux[,,ixflux]*1e6
ocnout=sum(ocntmp,na.rm=TRUE)
fosout=sum(fostmp,na.rm=TRUE)
firout=sum(firtmp,na.rm=TRUE)
#remove(list=c("ocnflux","fosflux","firflux","ocntmp","fostmp","firtmp","foot","dateall2"))
gc()
#--------------------------------------------------------------
# 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)
{
#if(lonarr[p]>=lonmin && lonarr[p]<=lonmax && latarr[p]>=latmin && latarr[p]<=latmax)
#{
#dy=as.integer(daylist[p])-dmin+1 #bug, who know if the dates are continuous or not
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
}
#i=ceiling((lonarr[p]-lonmin)/1.0+0.5) #check if the position is right
i=ceiling((lonarr[p]-lonmin)/3.0)
#nval=nval+1 #notice!!! this (i,j) is for NA
#notice: we maybe should use the data besides North America domain, not only nam1x1, but glb3x2...
#j=ceiling((latarr[p]-latmin)/1.0+0.5)
j=ceiling((latarr[p]-latmin)/2.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
}
# 3-hourly matching /agt/alt
hpos=ceiling((as.integer(hrlist[p])+0.000001+as.integer(milist[p])/60.0 )/3.0)
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
}
#}
}
#print(nval)
fbou=pbou/nval
#deltaco2=(recoout-gppout)+ocnout+fosout+firout #contributions from the biosphere and background fluxes
#fsimu=fbou+deltaco2
#################################################
# final results and output to files
#################################################
fnameout1=paste(outdir,"/",tb2,"/","samples_simulated.",tb2,"00","_",tb3,"00",".txt",sep="")
fn=fnameout1
outlist=NULL #c(inityr, initmo, initdy, inithr, initmi, inihgt, fbou)
for(i in 1:nsam)
{
deltaco2=neeoutarr[i]+ocnout+fosout+firout
fsimuarr[i]=fbou+deltaco2
outlist<-c(outlist, fsimuarr[i])
}
out1<-c(ident,eventid,round(outlist,4))
if(newfile)
write.table(t(out1), file=fnameout1, append=F, col.names=F, row.names=F, quote=F)
if(!newfile)
write.table(t(out1), file=fnameout1, append=T, col.names=F, row.names=F, quote=F)
newfile=F
}# end if foot NA
}#end loop mm
#-----------------------------------------------------------------------
# ouput rsults as *.nc files
#-----------------------------------------------------------------------
#source("/Users/weihe/Documents/Concensimulation/output.ncdf4.r")
#fn= "/Users/weihe/Documents/Concensimulation/resultstest/stilt_sib_ct_co2_2010-01-19_2010-01-29_simulated.txt"
fin=paste(fn,sep="") #already include eventid
fout=paste(substring(fn,1,nchar(fn)-4),".nc",sep="")
data=as.matrix(read.table(fin,header=FALSE))
nobs=dim(data)[1]
nmem=dim(data)[2]-2
vals=data[,2:(nsam+2)] #151 colums
write.results.netcdf(vals,nobs,nmem,fout) #write simulated results once
#-----------------------------------------------------------------------
#create stilt.ok file
file.create("stilt.ok")
# Time marker
TEnd=proc.time()
Tot<-TEnd-TBegin
print(Tot)
template.jb → templates/template.jb
+ 0
− 0
View file @ 9beec836
File moved
template.py → templates/template.py
+ 1
− 1
View file @ 9beec836
......@@ -27,7 +27,7 @@ sys.path.append(os.getcwd())
from da.tools.initexit import start_logger, validate_opts_args, parse_options, CycleControl
from da.tools.pipeline import ensemble_smoother_pipeline, header, footer
from da.platform.maunaloa import MaunaloaPlatform
from da.platform.capegrim import CapeGrimPlatform
from da.baseclasses.dasystem import DaSystem
from da.baseclasses.statevector import StateVector
from da.carbondioxide.obspack_globalviewplus2 import ObsPackObservations
......
template.rc → templates/template.rc
+ 0
− 0
View file @ 9beec836
File moved
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment