diff --git a/stilt.nc4.sib.co2.Feb.2015.May13.r b/stilt.nc4.sib.co2.Feb.2015.May13.r deleted file mode 100755 index 67b8af74100aa8ce25b9226f1b3b07aaeaf9599b..0000000000000000000000000000000000000000 --- a/stilt.nc4.sib.co2.Feb.2015.May13.r +++ /dev/null @@ -1,600 +0,0 @@ -# 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) diff --git a/template.jb b/templates/template.jb similarity index 100% rename from template.jb rename to templates/template.jb diff --git a/template.py b/templates/template.py similarity index 98% rename from template.py rename to templates/template.py index e672f37436a67db6ee66fc2569ebfaaace32d5a9..fd1ca9fba88395711bc2eae1227751a9a759b468 100755 --- a/template.py +++ b/templates/template.py @@ -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 diff --git a/template.rc b/templates/template.rc similarity index 100% rename from template.rc rename to templates/template.rc