R code about “ReddyPro” package
#The script was used to process original data from software "ReddyPro" by the package "REddyProc" in R.

#package must be installed by function "install.package("name of package")" if the package was not existed in R   #install.package("REddyProc")
#import required packages 
library(REddyProc)
library(segmented)
library(dplyr)
library(tidyr)
library(data.table)
library(lubridate)
library(timetk)

#reading data to process,which must have the customary name and form of variables and its unit in the following: "TIMESTAMP","NEE"(μmol/m2/s1),"H"(W/m2),"LE"(W/m2),"Ustar"(m/s),"VPD"(hPa),"rH"(%),"Tair"(℃),"Rg"(W/m2)
data <- read.csv("E:/arctic_research/analysis_sites/data_original/fc_CA-ARF.csv",header=TRUE,sep=',',stringsAsFactors = FALSE,na.strings=c("-9999","-9999.0"))
# data$TIMESTAMP <- parse_date_time(data$TIMESTAMP, "YmdHM") #recognize characters sequentially as year,month,day,hour,minute
data$TIMESTAMP <- as.POSIXct(data$TIMESTAMP,tz='UTC') #set it to time format
data <-  distinct(data,TIMESTAMP,.keep_all = TRUE) #data duplicated by specific column

#generate a continuous time series to make the data continuous in that the data processed by package "REddyProc" must be continuous 
cm_time <- as.integer((as.Date(range(data$TIMESTAMP,na.rm=TRUE)[2])-as.Date(range(data$TIMESTAMP,na.rm=TRUE)[1])+1))*48
cm_time <- strptime(paste(as.Date(range(data$TIMESTAMP,na.rm=TRUE)[1]),"00:00:00"),"%Y-%m-%d %H:%M:%S",tz="UTC")+30*60*1:cm_time
time <- as.data.frame(cm_time)
names(time)[1] <- "TIMESTAMP" 
data <- merge(time,data,by=c("TIMESTAMP"),all.x=TRUE) #merge 2 data in the same several or one columns 
data <- arrange(data,TIMESTAMP) #sort the data by the specific column

#divide day and nigh according to "Rg" or "PAR" or "H" and soon on; its demarcation point Rg > 20 W/m2 or PAR > 5 umol/m2/s1
if("dn" %in% names(data)){
  da <- data%>% subset(select=c(TIMESTAMP,NEE,dn))%>% arrange(TIMESTAMP)
  da$wi <- NA
  da$di <- NA
}else{
  if("Rg" %in% names(data)){ 
    da <- data%>% subset(select=c(TIMESTAMP,NEE,Rg))%>% arrange(TIMESTAMP)
    da$dn <- NA
    da$wi <- NA
    da$di <- NA
    da <- within(da,{
      dn[Rg > 10] <- 1
      dn[Rg <= 10] <- 0
    })
  }else{
    da <- data%>% subset(select=c(TIMESTAMP,NEE,H))%>% arrange(TIMESTAMP)
    da$dn <- NA
    da$wi <- NA
    da$di <- NA
    da <- within(da,{
      dn[H > -4.3] <- 1
      dn[H <= -4.3] <- 0
    })
  }
}
for(i in 1:nrow(da)){
  if(is.na(da[i,"dn"])){
    for(j in 1:round(nrow(da)/48)){
      da[i,"dn"] <- da[i+j*48,"dn"]
      if(is.na(da[i,"dn"])==F){
        break
      }
    }
    
  }
}#gap-filling day or nigh  used by later values at the same time
for(i in nrow(da):1){
  if(is.na(da[i,"dn"])){
    for(j in 1:round(nrow(da)/48)){
      da[i,"dn"] <- da[i-j*48,"dn"]
      if(is.na(da[i,"dn"])==F){
        break
      }
    }
    
  }
}#gap-filling day or nigh  used by earlier values at the same time

#eliminate outliers of NEE through MAD(median absolute deviation about the median)
for(i in 1:nrow(da)){
  da[i,"wi"] <- ceiling(i/(13*48)) #calculate window period at 13-day intervals
  if(i==1 | i==nrow(da)){
    da[i,"di"] <- NA
  }else{
    da[i,"di"] <- 2*da[i,"NEE"]-da[i-1,"NEE"]-da[i+1,"NEE"]
  }
}
d <- aggregate(da$di, by=list(dn=da$dn,wi=da$wi), FUN=median, na.rm=TRUE) #calculate the median of "di" at every day or night at every window period 
names(d)[3] <-c("md") 
da <- da %>% merge(d,by=c('dn','wi'),all.x=TRUE) %>% arrange(TIMESTAMP)
da$dmd <- abs(da$di-da$md) 
d <- aggregate(da$dmd, by=list(dn=da$dn,wi=da$wi), FUN=median, na.rm=TRUE);names(d) <-c("dn","wi","mad") 
da <- da %>% merge(d,by=c('dn','wi'),all.x=TRUE) %>% arrange(TIMESTAMP)%>%
  mutate(up=md+((4*mad)/0.6745),low=md-((4*mad)/0.6745))%>% arrange(TIMESTAMP) #generate reasonable interval of NEE. constant: 4,5.5,7, the default is 4
da <- within(da,{
  NEE[di > up | di < low] <- NA
})#eliminate outliers of NEE according to lower and up threshold
data$NEE <- da$NEE

#calculate Ustar(friction wind speed at night)
data <- plyr::rename(data, c("TIMESTAMP"="DateTime")) #rename "TIMESTAMP" to "DateTime"
lat=33.2;lon=96.55;tz=7 #set latitude, longitude and time zone of corresponding site
edd.c <- sEddyProc$new('SERVIRST',data,ColPOSIXTime="DateTime",LatDeg=lat,LongDeg=lon,TimeZoneHour=tz,names(data[-1])) #set necessary scenarios to calculate Ustar
ustarth <- edd.c$sEstUstarThresholdDistribution(nSample=200L,probs=c(0.05,0.25,0.50,0.95)) #set sample size
ustarth%>%filter(aggregationMode=="season")%>%select(uStar,"5%","25%","50%","95%") #extract seasonal Ustar
ustarthseason <- usGetSeasonalSeasonUStarMap(ustarth)
ustarfix <- colnames(ustarthseason)[-1]
edd.c$sSetUstarScenarios(ustarthseason,uStarSuffixes=ustarfix) #set Ustar scenarios corresponding gap-filling NEE

#impute NEE
edd.c$sMDSGapFillUStarScens('NEE',FillAll=TRUE)

#impute each environment variable,except NEE, in the same method as NEE in order to split NEE into Re and GPP
edd.c$sMDSGapFill('Tair',FillAll=TRUE)
edd.c$sMDSGapFill('rH',FillAll=TRUE)
edd.c$sMDSGapFill('VPD',FillAll=TRUE)
edd.c$sMDSGapFill('Rg',FillAll=TRUE)
edd.c$sMDSGapFill('H',FillAll=TRUE)
edd.c$sMDSGapFill('LE',FillAll=TRUE)

#split NEE into Re and GPP
resP <- lapply(ustarfix,function(ustarfix){edd.c$sMRFluxPartition(suffix=ustarfix)})

# #paint result processed by "REddyProc"
# grep("NEE.*_f$|Reco|GPP",names(edd.c$sExportResults()),value=TRUE) #index variable according to specific requirements
# edd.c$sPlotDailySumsY("NEE_uStar_f",Year=2011)
# edd.c$sPlotFingerprintY("NEE_uStar_f",Year=2011)
# edd.c$sPlotHHFluxesY("NEE_uStar_f",Year=2011)

# extract and output final result according to your requirements
ld <- cbind(data,edd.c$sExportResults())  #merge original data and result from package "REddyProc" processed in row
ld <- ld[,c("DateTime","NEE_uStar_f","NEE_uStar_fqc")] #extract need variable according to your requirements
data <- plyr::rename(ld, c("DateTime"="TIMESTAMP","NEE_uStar_f"="NEE_f","NEE_uStar_fqc"="NEE_qc")) #rename variable
data$month <- month(data$TIMESTAMP) #generate month variable from TIMESTAMP
data$year <- year(data$TIMESTAMP) #generate year variable from TIMESTAMP
data <- within(data,{NEE_f[NEE_qc > 1] <- NaN}) #make NEE_f be the missing value on the condition that the quality of NEE is classify as 2(worst quality of NEE)
dm <- aggregate(data$NEE_f, by=list(year=data$year,month=data$month), FUN=mean, na.rm=TRUE) #calculate the annual and monthly average of NEE by grouping
dm <- plyr::rename(dm, c("x"="NEE_f"))
n<-paste('E:/Arctic_study/data/','name of data','.csv',sep='') #output path and name of data
write.csv(dm,file = n,row.names=F) #write data