############################################################
#CALCULATION OF THE INTEGRAL OVER THE UNIVERSAL FUNCTION
  #after Businger (1971) in the form of HC6gstrC6m (1988)
  #from Foken (2006) Eq. (2.85) - (2.88)
############################################################

univfun <- function(
  zmeas,
  Lvirt,
  sigma_range=c(-2, 1)
) {

  #calculation of local stability sigma [-]
  sigma <- zmeas / Lvirt
  
  #assign a flag for universal function being calculated within defined stability range
  sigma_flag <- 0
  
  #instable case
  if(sigma < 0) {
    
    #catch cases for which parametrization is not defined
    if (sigma < sigma_range[1]) { 
      sigma <- sigma_range[1]
      sigma_flag <- -1
    }
    
    #actual calculation
    x <- (1 - 19.3 * sigma)^(1/4)
    y <- 0.95 * (1 - 11.6 * sigma)^(1/2)
    Psi <- log(((1 + x^2) / 2) * ((1 + x) / 2)^2) - 2 * atan(x) + pi / 2
    rm(x, y)
    
    #stable case
  } else {
    
    #catch cases for which parametrization is not defined
    if (sigma > sigma_range[2]) {
      sigma <- sigma_range[2]
      sigma_flag <- 1
    }
    
    #actual calculation
    Psi <- -6 * sigma
    
  }
  
  #return result
  OUT <- list(Psi=Psi, sigma_flag=sigma_flag)
  return(OUT)
  
  #universal function for momentum exchange after Skeib (1980)
  #not integrated, only defined for sigma_range <- c(-2,2)
  #only defined for sigma_range <- c(-2,2)
  
  #any case not fulfilling sigma_range
  #Psi <- NA
  #labile case
  #if(sigma > sigma_range[1] & sigma < -0.0625) Psi <- (-sigma / 0.0625)^(-1/4)
  #neutral case
  #if(sigma >= -0.0625 & sigma <= 0.125) Psi <- 1
  #stable case
  #if(sigma > 0.125 & sigma < sigma_range[2]) Psi <- (sigma / 0.125)
  
}


############################################################
#CALCULATION OF AERODYNAMIC ROUGHNESS LENGTH z0 (NEEDED FOR K04 footprint model)
############################################################

roughness <- function(
  zmeas,
  Lvirt,
  hor,
  ustar,
  zLrange=c(-2, 1)
) {

  dum <- sapply(1:length(Lvirt), function(x) 
    #x <-1
    univfun(zmeas=zmeas[x],Lvirt=Lvirt[x],sigma_range=zLrange
  )
  )
  dum <- data.frame(t(dum))
  Psi <- unlist(dum$Psi)
  sigma_flag <- unlist(dum$sigma_flag)
  rm(dum)
  
  #calculation of roughness length d_z_0 [m]
  d_z_0 <- zmeas / exp(kap * hor / ustar + Psi)
  names(d_z_0) <- "d_z_0"
  
  #return result
  return(d_z_0)

}



############################################################
#STATIONARITY TESTS
############################################################

REYNstat <- function(
  data,			#data frame with EC data
  type=c(1, 2, 3)[2],	#analysis with trend (1) or internal stationarity (2) or both (3)
  whr_flux, #for which fluxes to perform stationarity test
  NOsusa=6,		#number of subsamples for trend=F
  FcorMCO,
  FcorPOT,
  FcorPOTl=FcorPOTl,
  FcorWPL=FcoWPL
) {


#-----------------------------------------------------------
#BASIC SETUP

  #fluxes including trend
    tren <- REYNflux(
      data=data,
      basetype="mean",
      FcorMCO=FcorMCO,
      FcorPOT=FcorPOT,
      FcorWPL=FcorWPL
    )

###
if(type %in% c(1, 3)) {
###

#-----------------------------------------------------------
#TREND EFFECT

  #fluxes after trend removal
    detr <- REYNflux(
      data=data,
      basetype="trend",
      FcorMCO=FcorMCO,
      FcorPOT=FcorPOT,
      FcorWPL=FcorWPL,
      FcorPOTl=FcorPOTl
    )

  #deviation [%]
    crit_t <- ((detr$mn - tren$mn) / tren$mn * 100)[whr_flux]

  #clean up
    rm(detr)


###
} else crit_t <- NULL
###


###
if(type %in% c(2, 3)) {
###

#-----------------------------------------------------------
#INTERNAL INSTATIONARITIES

  #class boundaries
    sampBO <- round(seq(1, nrow(tren$data), length.out=NOsusa + 1))
    sampBO[length(sampBO)] <- sampBO[length(sampBO)] + 1

  #list with indexes of subsamples
    sampLU <- sapply(1:(length(sampBO) - 1), function(x) seq(sampBO[x], sampBO[x + 1] - 1))

  #results for the subsamples
    sampRE <- sapply(1:NOsusa, function(x) 
      #x <- 6
      REYNflux(
        data=data[sampLU[[x]],],
	      basetype="mean",
	      FcorMCO=FcorMCO,
	      FcorPOT=FcorPOT,
	      FcorPOTl=FcorPOTl,
	      FcorWPL=FcorWPL
      )$mn[,whr_flux]
    )
    sampRE <- data.frame(matrix(unlist(sampRE), ncol=length(whr_flux), byrow=T))
    dimnames(sampRE)[[2]] <- whr_flux

  #stationarity criteria
    crit_i <- (colMeans(sampRE) - tren$mn[whr_flux]) / tren$mn[whr_flux] * 100

  #clean up
    rm(tren, NOsusa, sampBO, sampLU, sampRE)

###
} else crit_i <- NULL
###


#-----------------------------------------------------------
#AGGREGATE AND RETURN RESULTS

  #aggregate results
    crit <- list()
    if(!is.null(crit_t)) crit$trend=crit_t
    if(!is.null(crit_i)) crit$subsa=crit_i

  #return results
    return(crit)

}



############################################################
#INTEGRAL TURBULENCE CHARACTERISTICS
############################################################


REYNitcs <- function(
  sigma,
  latitude,
  stdev,
  scale
) {

  
  #-----------------------------------------------------------
  #MODELLED ITCS
  
  #STANDARD MODEL
  #assign model coefficients
  #stdevU / u_star
  ITCuC <- if(sigma < -0.032) c(4.15, 1/8) else c(2.7, 0)
  #stdevW / u_star
  ITCwC <- if(sigma < -0.032) c(2.0, 1/8) else c(1.3, 0)
  #stdevT / T_star_SL
  ITCtC <-  if(sigma < -1) c(1, -1/3) else 
    if(sigma >= -1 & sigma < -0.062) c(1, -1/4) else
      if(sigma >= -0.062 & sigma < 0.02) c(0.5, -1/2) else
        if(sigma >= 0.02) c(1.4,  -1/4)
  #calculate models; the absolute value of stability is used, since negative values often result in NAs of the potency
  ITCuM <- ITCuC[1] * abs(sigma)^ITCuC[2]
  ITCwM <- ITCwC[1] * abs(sigma)^ITCwC[2]
  ITCtM <- ITCtC[1] * abs(sigma)^ITCtC[2]
  
  #UPDATED MODEL (THOMAS AND FOKEN, 2002)
  #coriolis parameter
  f <- coriolis(latitude)
  #update models
  if(sigma > -0.2 & sigma < 0.4) {
    ITCuM <- 0.44 * log(f / scale$ustar) + 6.3
    ITCwM <- 0.21 * log(f / scale$ustar) + 3.1
  }
  #clean up
  rm(ITCuC, ITCwC, ITCtC, f)
  
  
  #-----------------------------------------------------------
  #COMPARE TO MEASURED ITCS
  
  #calculate observed ITCs
  ITCuO <- stdev$hor_x / scale$ustar
  ITCwO <- stdev$hor_z / scale$ustar
  ITCtO <- stdev$T_air / scale$Tstar
  
  #final criteria [%]
  #individual
  ITC <- data.frame(
    hor_x=(abs(ITCuO - ITCuM) / ITCuM) * 100,
    hor_z=(abs(ITCwO - ITCwM) / ITCwM) * 100,
    T_air=(abs(ITCtO - ITCtM) / ITCtM) * 100
  )
  
  #combined
  ITC$ustar=sqrt(ITC$hor_x^2 + ITC$hor_z^2)
  ITC$F_H_kin=sqrt(ITC$hor_z^2 + ITC$T_air^2)
  
  #clean up
  rm(ITCuO, ITCuM, ITCwO, ITCwM, ITCtO, ITCtM)
  
  #-----------------------------------------------------------
  #EXPORT RESULTS
  
  return(ITC)
}



############################################################
#INTEGERAL LENGTH SCALES
############################################################
#The turbulence length scale, l , is a physical quantity describing the size of the large energy-containing eddies in a turbulent flow.

INTsca <- function(
  parcel,
  depe
) {

  #fill gaps via linear interpolation
    depe <- approx(parcel, depe, xout=parcel)[[2]]

  #get rid of NAs at start and end
    dum_NA <- na.omit(data.frame(parcel=parcel, depe=depe))
    parcel <- dum_NA$parcel
    depe <- dum_NA$depe
    rm(dum_NA)

  #demeaning and detrending
    depe <- lm(depe ~ parcel)$residuals

  #path through air [m] per increment, the stepwidth of the integral scale        #del:
  incr <- max(parcel - min(parcel)) / length(parcel)
  
  #calculate auto-correlation function
    lag <- 10; crit <- 1
    while(crit > 0) {
      lag <- lag * 2
      ACF <- acf(depe, lag.max = lag, type = "correlation", plot = FALSE, na.action = na.fail, demean = TRUE)
      crit <- min(ACF$acf)
    }
    
    require(EMD)

  #integral length scale: typical size of the largest or most energy-transporting eddies.
  #a) integral over acf until first zero crossing (Bange, 2002);
  #b) first maximum of the integral (Lenschow, 1986);
    #find first zero crossing
#      whr_zero <- ACF$lag[zeroCross(ACF$acf, slope="negative")[1]]
      whr_zero <- extrema(y=c(ACF$acf, sin(1:10)))$cross[1,2]
    
    
    #for each cell, weight distance increment with correlation coefficient
      I <- sum(ACF$acf[1:whr_zero]) * incr
    #alternatively: all in one, assume that correlation monotonously decreases after peak
      #I <- max(cumsum(ACF$acf)) * incr

  #return result
    return(I)

}



############################################################
#STATISTICAL ERROR FOR FLUXES
############################################################

REYNerro <- function(
  isca,		#integral scale lengths
  mn,		#mean values
  rc,		#flux correlation coefficient
  Lf		#averaging distance (flight length)
) {

  
  #-----------------------------------------------------------
  #VARIABLES
  
  #scalar length scales
  ls <- isca$scal
  
  #lookup table for matching scalar length scales to flux correlation coefficients
  whr_ls <- rbind(  #change!
    c("ustar2_x", "hor_x"),
    c("ustar2_y", "hor_y"),
    c("H", "fst_Tair"), 
    c("E", "fst_FD_mole_H2O_insitu"),  #I used FD_mole_H2O_flux, you used E_en, does it matter?
    c("out_CH4_flux",  "FD_mole_CH4_con"),  #I used FD_mole_CH4_flux, you used F_CH4_mass, does it matter?
    c("out_CO2_flux","fst_FD_mole_CO2_insitu")
  )
  
#   str(ls)
#   'data.frame':  1 obs. of  8 variables:
#     $ hor_x                 : num 84.8
#   $ hor_y                 : num 186
#   $ hor_z                 : num 57.4
#   $ fst_Tair              : num 56.7
#   $ fst_SONIC_T           : num 76.4
#   $ fst_FD_mole_H2O_insitu: num 115
#   $ fst_FD_mole_CO2_insitu: num 1.6
#   $ FD_mole_CH4_con       : num 23.5
#   
  #variance length scales
  lv <- isca$vari
  #flux length scales
  lf <- isca$flux
  #matching flux correlation coefficients
  rc <- rc[match(dimnames(lf)[[2]], dimnames(rc)[[2]])]
  
  
  #-----------------------------------------------------------
  #RANDOM ERROR
  
  #variance (Lenschow, 1994 Eq. 36)
  erra_vari <- sqrt(2 * lv / Lf) * 100
  
  #fluxes
  
  #directly from flux length scales (Lenschow, 1986 Eq. (7) == Lenschow, 1994 Eq. (48))
  #error estimate
  #large error in F_CH4_mass due to four times lower correlation coefficient
  erra_flux <- data.frame(sqrt(2 * lf / Lf * (rc^(-2) + 1)) * 100)
  #error reproduction
  #Bange: u_star2_x and u_star2_y are orthogonal == indepenend -> for random error we can use Gauss reproduction
  erra_flux$ustar <- (
    (mn$ustar2_x * erra_flux$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
      (mn$ustar2_y * erra_flux$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
  )^(1/4) / mn$ustar * 100
  
  
  
  #upper limit from scalar length scales
  #upper limt for lenght scale following Mann (1994) Eq. (9), Bange (2002) Eq. (10)
  
  # situation that dimnames(rc)[[2]] has "ustar" and "FD_mole_H2O_flux" while whr_ls doesn't. 
  if (dimnames(rc)[[2]][3] == "ustar") {
    rc <- rc[,-3]
  }
    
  if (dimnames(rc)[[2]][4] == "E") {
    rc <- rc[,-4]
  } 
  
  lf_max <- data.frame(sapply(1:length(rc), function(x) {
    #get corresponding scalar length scale from conversion table
    #x <- 6
   if (is.na(match(dimnames(rc)[[2]][x],whr_ls[,1])) != TRUE  ) {
     if(is.na(match(whr_ls[match(dimnames(rc)[[2]][x],whr_ls[,1]),2],dimnames(ls)[[2]])) != TRUE){
          ls_loc <- ls[,whr_ls[match(dimnames(rc)[[2]][x],whr_ls[,1]),2]]
          sqrt(ls$hor_z * ls_loc) / abs(rc[x])
     }
   }
    else{
      ls_loc <- NaN
      NaN}
    
  }
  ))
  #error estimate
  #with Mann / Bange length scale
  errm_flux <- data.frame(sapply(1:length(rc), function(x)  sqrt(2 * lf_max[x] / Lf * (rc[x]^(-2) + 1)) * 100))
    
  #according to Lenschow, 1994 Eq. (49); Bange (2002) Eq. (16) says: WRONG
  #errm <- data.frame(sapply(1:length(rc), function(x) abs(2 / rc[x] * sqrt(min(ls$hor_z, ls[,x]) / Lf)) * 100))
  #error reproduction (Bange, 2007 Eq.(3.32))
  #u_star: Bange u_star2_x and u_star2_y are orthogonal == indepenend -> for random error we can use Gauss reproduction
  errm_flux$ustar <- (
    (mn$ustar2_x * errm_flux$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
      (mn$ustar2_y * errm_flux$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
  )^(1/4) / mn$ustar * 100
  
  #save to list
  ran <- list()
  #sampling error in variance
  ran$vari$act=erra_vari
  #actual flux random error (flux length scales)
  ran$flux$act=erra_flux		
  #maximum flux random error (scalar length scales)
  ran$flux$max=errm_flux		
  
  #clean up
  rm(erra_flux, erra_vari, errm_flux)
  
  
  #-----------------------------------------------------------
  #SYSTEMATIC ERROR
  
  #variance (Lenschow, 1994 Eq. 36)
  ersa_vari <- 2 * lv / Lf * 100
  
  #fluxes
  
  #directly from flux length scales (Lenschow, 1994 Eq. (27) == Bange, 2007 Eq. (3.20))
  #error estimate
  ersa_flux <- data.frame(2 * lf / Lf * 100)
  #error reproduction
  #u_star upper bound (no independence and randomness required)
  ersa_flux$ustar <- (
    (mn$ustar2_x * ersa_flux$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
      (mn$ustar2_y * ersa_flux$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
  )^(1/4) / mn$ustar * 100
  
  #upper limit from scalar length scales
  #error estimate
  #with Mann / Bange length scale
#? THis is length scale, why it is ersm_flux?
  ersm_flux <- data.frame(sapply(1:length(rc), function(x) 2 * lf_max[x] / Lf * 100))
  #according to Lenschow, 1994 Eq. (29); same!
  #ersm <- data.frame(sapply(1:length(rc), function(x) abs(2 / rc[x] * sqrt(ls$w_hor * ls[,x] ) / Lf) * 100))
  #error reproduction
  #u_star upper bound (no independence and randomness required)
  ersm_flux$ustar <- (
    (mn$ustar2_x * ersm_flux$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
      (mn$ustar2_y * ersm_flux$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
  )^(1/4) / mn$ustar * 100
  
  #save to list
  sys <- list()
  #sampling error in variance
  sys$vari$act=ersa_vari
  #actual flux systematic error (flux length scales)
  sys$flux$act=ersa_flux
  #maximum flux systematic error (scalar length scales)
  sys$flux$max=ersm_flux
  
  #clean up
  rm(ls, lv, lf, ersa_flux, ersa_vari, ersm_flux, lf_max)
  
  
  #-----------------------------------------------------------
  #AGGREGATE AND EXPORT RESULTS
  
  #aggregate
  export <- list(
    ran=ran,
    sys=sys
  )
  #clean up
  rm(ran, sys)
  #export
  return(export)

}



############################################################
#COMPLETE FLUX COMPUTATION SEQUENCE
############################################################

REYNcomp <- function(
  eddy.data=eddy.data,
  whr_scal=c("hor_x", "hor_y", "hor_z", "fst_Tair", "fst_SONIC_T", "fst_FD_mole_H2O_insitu", "fst_FD_mole_CO2_insitu","FD_mole_CH4_con"),  #change
  whr_flux=c("ustar2_x", "ustar2_y", "ustar", "H", "E", "FD_mole_CH4_flux","FD_mole_H2O_flux","FD_mole_CO2_flux"), #change these variables' names
  latitude=tower_latitude,
  basetype=c("mean", "trend", "3order")[1],
  stattype=c(1, 2, 3)[1],
  NOsusa=6,
  FcorMCO=TRUE,  #change
  FcorPOT=TRUE,
  FcorPOTl=NULL,
  FcorWPL=FALSE #change
) {


#-----------------------------------------------------------
#CALCULATE FLUXES
  REYN <- REYNflux(
    data=eddy.data,
    basetype=basetype,
    FcorMCO=FcorMCO,
    FcorPOT=FcorPOT,
    FcorPOTl=FcorPOTl,
    FcorWPL=FcorWPL
  )


#-----------------------------------------------------------
#CALCULATE AERODYNAMIC ROUGHNESS LENGTH
  REYN$mn$z0 <- roughness(
    zmeas=REYN$mn$uls_z,
    Lvirt=REYN$mn$Lvirt,
    hor=REYN$mn$hor_x,
    ustar=REYN$mn$ustar
  )


#-----------------------------------------------------------
#STATIONARITY CRITERIA [%]
  REYN$stat <- REYNstat(
    data=eddy.data,
    type=c(1, 2, 3)[stattype],
    whr_flux=whr_flux,
    NOsusa=NOsusa,
    FcorMCO=FcorMCO,
    FcorPOT=FcorPOT,
    FcorPOTl=FcorPOTl,
    FcorWPL=FcorWPL
  )


#-----------------------------------------------------------
#INTEGRAL TURBULENCE CHARACTERISTICS [%]
  #data frame for deviations
    stdev <- data.frame(
      hor_x=REYN$sd$hor_x,
      hor_z=REYN$sd$hor_z,
      T_air=REYN$sd$fst_Tair
    )

  #data frame for scales
    scale <- data.frame(
      ustar=REYN$mn$ustar,
      Tstar=REYN$mn$Tstar
    )

  #calculation
    REYN$itcs <- REYNitcs(
      sigma=REYN$mn$zLvirt,   #stability
      latitude=latitude,
      stdev=stdev,
      scale=scale
    )

  #clean up
    rm(stdev, scale)


#-----------------------------------------------------------
#INTEGRAL LENGTH SCALES
  #available variables in REYN$imfl
#scalar length scales
#calculate scales

isca_scal <- sapply(whr_scal, function(x) 
  {if(length(which(is.na(REYN$imfl[,x]))) != length(REYN$imfl[,x])){
    INTsca(parcel=REYN$data$parcel, depe=REYN$imfl[,x])
  }else {
    NaN
  }
}
)

isca_scal <- as.data.frame( matrix(isca_scal, ncol=length(whr_scal)) )
attributes(isca_scal)$names <- whr_scal

#variance length scales
#calculate scales
isca_vari <- sapply(whr_scal, function(x) 
  {if(length(which(is.na(REYN$imfl[,x]))) != length(REYN$imfl[,x])){
      INTsca(
      parcel=REYN$data$parcel,
      depe=REYN$imfl[,x]^2
      )
      }else {
        NaN
      }
  }
)
isca_vari <- as.data.frame( matrix(isca_vari, ncol=length(whr_scal)) )
attributes(isca_vari)$names <- whr_scal

#flux length scales
#exclude data with > 10% immidiate flucuations not available, e.g. u_star;
#sampling error in u_star has to be assessed by its both components u_star_x and u_star_y
whr_not <- sapply(whr_flux, function(x) length(which(is.na(REYN$imfl[,x]))) / nrow(REYN$imfl))
whr_not <- which(whr_not > 0.5)
if(length(whr_not) == 0) whr_flux_loc <-  whr_flux else whr_flux_loc <-  whr_flux[-whr_not]

#calculate scales    
isca_flux <- sapply(whr_flux_loc, function(x) 
 # x <- "out_CH4_flux"
  INTsca(
  parcel=REYN$data$parcel,
  depe=REYN$imfl[,x]
)
)
isca_flux <- as.data.frame( matrix(isca_flux, ncol=length(whr_flux_loc)) )
attributes(isca_flux)$names <- whr_flux_loc

#assign to list
REYN$isca <- list(
  scal=isca_scal,
  vari=isca_vari,
  flux=isca_flux
)

#clean up
rm(isca_scal, isca_vari, isca_flux, whr_flux_loc, whr_not)

#-----------------------------------------------------------
#STATISTICAL ERROR FOR FLUXES
  REYN$erro <- REYNerro(
    isca=REYN$isca,		#integral scale lengths
    mn=REYN$mn,			#mean values
    rc=REYN$cor[,-which(dimnames(REYN$cor)[[2]]=="Ht")],		#flux correlation coefficient
    Lf=REYN$max$parcel
  )


#-----------------------------------------------------------
#RETURN RESULT
  return(REYN)

}