############################################################
#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
) {

#calculation of local stability zeta [-]
  zeta <- zmeas / Lvirt

#instable case
  if(zeta < 0) {
    
    #catch cases for which parametrization is not defined
      if (zeta < -2) zeta <- -2
      
    x <- (1 - 19.3 * zeta)^(1/4)
    y <- 0.95 * (1 - 11.6 * zeta)^(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 (zeta > 1) zeta <- 1

    Psi <- -6 * zeta
  }

#return result
  return(Psi)

#universal function for momentum exchange after Skeib (1980)
  #not integrated, only defined for zLrange <- c(-2,2)
  #only defined for zLrange <- c(-2,2)

  #any case not fulfilling zLrange
    #Psi <- NA
  #labile case
    #if(zeta > zLrange[1] & zeta < -0.0625) Psi <- (-zeta / 0.0625)^(-1/4)
  #neutral case
    #if(zeta >= -0.0625 & zeta <= 0.125) Psi <- 1
  #stable case
    #if(zeta > 0.125 & zeta < zLrange[2]) Psi <- (zeta / 0.125)

}


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

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

#von Karman constant
  kappa <- 0.4

#calculation of local stability zeta [-]
  zeta <- zmeas / Lvirt

#integral over the universal function after Businger (1971) in the form of HC6gstrC6m (1988)
#Foken (2006) Eq. (2.85) - (2.88)
  Psi <- sapply(1:length(Lvirt), function(x) univfun(zmeas[x], Lvirt[x]) )

#calculation of roughness length z0 [m]
  z0 <- zmeas / exp(kappa * hor / ustar + Psi)

#return result
  return(z0)

}



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

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


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

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

   vars <- attributes(tren$mn)$names
  #which fluxes to compare
    whr_flux <- sapply(1:length(flux_vars), function(x) which(vars == flux_vars[x]))


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

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

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

  #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) REYNflux(
      #x <- 1
	data=tren$mcda[sampLU[[x]],],
	basetype="mean",
	FcorMCO=F,
	FcorPOT=FcorPOT,
	FcorPOTl=tren$mn$psta,
	FcorWPL=FcorWPL
      )$mn[,whr_flux]
    )
    sampRE <- data.frame(matrix(unlist(sampRE), ncol=length(whr_flux), byrow=T))
    attributes(sampRE)$names <- attributes(tren$mn)$names[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(
  zLvirt,
  latitude,
  sigma,
  scale
) {


#-----------------------------------------------------------
#MODELLED ITCS

  #STANDARD MODEL

    #assign model coefficients
      #sigmaU / ustar
	ITCuC <- if(zLvirt < -0.032) c(4.15, 1/8) else c(2.7, 0)

      #sigmaW / ustar
	ITCwC <- if(zLvirt < -0.032) c(2.0, 1/8) else c(1.3, 0)

      #sigmaT / Tstar
	ITCtC <-  if(zLvirt < -1) c(1, -1/3) else 
		  if(zLvirt >= -1 & zLvirt < -0.062) c(1, -1/4) else
		  if(zLvirt >= -0.062 & zLvirt < 0.02) c(0.5, -1/2) else
		  if(zLvirt >= 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(zLvirt)^ITCuC[2]
      ITCwM <- ITCwC[1] * abs(zLvirt)^ITCwC[2]
      ITCtM <- ITCtC[1] * abs(zLvirt)^ITCtC[2]


  #UPDATED MODEL (THOMAS AND FOKEN, 2002)

    #coriolis parameter
      f <- coriolis(latitude)

    #update models
      if(zLvirt > -0.2 & zLvirt < 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 <- sigma$u / scale$ustar
    ITCwO <- sigma$w / scale$ustar
    ITCtO <- sigma$Temp / scale$Tstar

  #final criteria [%]
    #individual
      ITC <- data.frame(
	u=(abs(ITCuO - ITCuM) / ITCuM) * 100,
	w=(abs(ITCwO - ITCwM) / ITCwM) * 100,
	Temp=(abs(ITCtO - ITCtM) / ITCtM) * 100
      )

    #combined
      ITC$ustar=sqrt(ITC$u^2 + ITC$w^2)
      ITC$H=sqrt(ITC$w^2 + ITC$Temp^2)
      ITC$E=ITC$w
      ITC$C=ITC$w

  #clean up
    rm(ITCuO, ITCuM, ITCwO, ITCwM, ITCtO, ITCtM)

#-----------------------------------------------------------
#EXPORT RESULTS

  return(ITC)

}



############################################################
#INTEGERAL LENGTH SCALES
############################################################


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 <- 1e2; crit <- 1
    while(crit > 0) {
      lag <- lag * 2
      ACF <- acf(depe, lag.max = lag, type = "correlation", plot = F, na.action = na.fail, demean = TRUE)
      crit <- min(ACF$acf)
    }

  #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 <- ACF$lag[nearest(x=ACF$acf, xval=0)[1]]
    
    
    #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)
  wplc,
  FcorWPL
) {


#-----------------------------------------------------------
#VARIABLES

  #available variables in isca
    vars <- attributes(isca)$names

  #scalar length scales !!must match the order of the fluxes!!
    scal <- c("hor_x", "hor_y", "fst_SONIC_T", "fst_FD_mole_H2O_insitu", "fst_FD_mole_CO2_insitu","FD_mole_CH4_con")
    whr_scal <- sapply(1:length(scal), function(x) which(vars == scal[x]))
    ls <- isca[,whr_scal]

  #variance length scales
    whr_vari <- which(vars == "hor_x") : which(vars == "FD_mole_CH4_con_var")
    lv <- isca[,whr_vari]

  #flux length scales
    whr_flux <- which(vars == "ustar2_x") : which(vars == "E")
    lf <- isca[,whr_flux]

  #absolute weights in Gaussian error reproduction
    if(FcorWPL == T) wplc <- abs(wplc)

  #clean up	
    rm(vars, scal, whr_scal, whr_vari, whr_flux)


#-----------------------------------------------------------
#RANDOM ERROR

  #variance (Lenschow, 1994 Eq. 36)
    erra <- sqrt(2 * lv / Lf) * 100

  #fluxes

    #directly from flux length scales (Lenschow, 1986 Eq. (7) == Lenschow, 1994 Eq. (48))

      #error estimate
	erra <- data.frame(c(erra, sqrt(2 * lf / Lf * (rc^(-2) + 1)) * 100))

      #error reproduction
	#Bange: ustar2_x and ustar2_y are orthogonal == indepenend -> for random error we can use Gauss reproduction
	  erra$ustar <- (
	    (mn$ustar2_x * erra$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
		(mn$ustar2_y * erra$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
	  )^(1/4) / mn$ustar * 100

	#for WPL corrected fluxes (upper bound)
	  if(FcorWPL == T) {
	      erra$E <- abs(wplc$H2E * erra$H) + abs(wplc$E2E * erra$E)
	      erra$C <- abs(wplc$H2C * erra$H) + abs(wplc$E2C * erra$E) + abs(wplc$C2C * erra$C)
	  }

    #upper limit from scalar length scales

      #upper limt for lenght scale following Mann (1994) Eq. (9), Bange (2002) Eq. (10)
	lf_max <- data.frame(sapply(1:length(rc), function(x) sqrt(isca$hor_z * ls[,x]) / abs(rc[x])))

      #error estimate
	#with Mann / Bange length scale
	  errm <- 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(isca$hor_z, ls[,x]) / Lf)) * 100))

      #error reproduction (Bange, 2007 Eq.(3.32))
	#ustar: Bange ustar2_x and ustar2_y are orthogonal == indepenend -> for random error we can use Gauss reproduction
	  errm$ustar <- (
	    (mn$ustar2_x * errm$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
		(mn$ustar2_y * errm$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
	  )^(1/4) / mn$ustar * 100

	#for WPL corrected fluxes (upper bound)
	  if(FcorWPL == T) {
	      errm$E <- abs(wplc$H2E * errm$H) + abs(wplc$E2E * errm$E)
	      errm$C <- abs(wplc$H2C * errm$H) + abs(wplc$E2C * errm$E) + abs(wplc$C2C * errm$C)
	  }

    #save to dataframe
      ran <- list(
	act=erra,		#actual random error (flux length scales)
	max=errm		#maximum random error (scalar length scales)
      )

    #clean up
      rm(erra, errm)


#-----------------------------------------------------------
#SYSTEMATIC ERROR

  #variance (Lenschow, 1994 Eq. 36)
    ersa <- 2 * lv / Lf * 100

  #fluxes

    #directly from flux length scales (Lenschow, 1994 Eq. (27) == Bange, 2007 Eq. (3.20))

      #error estimate
	ersa <- data.frame(c(ersa, 2 * lf / Lf * 100))

      #error reproduction

	#ustar upper bound (no independence and randomness required)
	  ersa$ustar <- (
	    (mn$ustar2_x * ersa$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
		(mn$ustar2_y * ersa$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
	  )^(1/4) / mn$ustar * 100

	#for WPL corrected fluxes (upper bound)
	  if(FcorWPL == T) {
	      ersa$E <- abs(wplc$H2E * ersa$H) + abs(wplc$E2E * ersa$E)
	      ersa$C <- abs(wplc$H2C * ersa$H) + abs(wplc$E2C * ersa$E) + abs(wplc$C2C * ersa$C)
	  }

    #upper limit from scalar length scales

      #error estimate

	#with Mann / Bange length scale
	  ersm <- 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(isca$hor_z * ls[,x] ) / Lf) * 100))

      #error reproduction

	#ustar upper bound (no independence and randomness required)
	  ersm$ustar <- (
	    (mn$ustar2_x * ersm$ustar2_x / 100 * mn$ustar2_x / mn$ustar^2)^2 +
		(mn$ustar2_y * ersm$ustar2_y / 100 * mn$ustar2_y / mn$ustar^2)^2
	  )^(1/4) / mn$ustar * 100

	#for WPL corrected fluxes (upper bound)
	  if(FcorWPL == T) {
	      ersm$E <- abs(wplc$H2E * ersm$H) + abs(wplc$E2E * ersm$E)
	      ersm$C <- abs(wplc$H2C * ersm$H) + abs(wplc$E2C * ersm$E) + abs(wplc$C2C * ersm$C)
	  }

    #save to dataframe
      sys <- list(
	act=ersa,		#actual systematic error (flux length scales)
	max=ersm		#maximum systematic error (scalar length scales)
      )

    #clean up
      rm(ls, lv, lf, ersa, ersm, 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,
  latitude=tower_latitude,
  basetype=basetype,
  stattype=c(1, 2, 3)[2],
  NOsusa=6,
  FcorMCO=FcorMCO,
  FcorPOT=FcorPOT,
  FcorPOTl=NULL,
  FcorWPL=FcorWPL
) {


#-----------------------------------------------------------
#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,
    zLrange=c(-Inf, Inf)
  )


#-----------------------------------------------------------
#STATIONARITY CRITERIA [%]
  REYN$stat <- REYNstat(
    data=eddy.data,
    flux_vars=c("Ht", "FD_mole_CH4_flux","FD_mole_CO2_flux","FD_mole_H2O_flux","E"),
    type=c(1, 2, 3)[stattype],
    NOsusa=NOsusa,
    FcorMCO=FcorMCO,
    FcorPOT=FcorPOT,
    FcorWPL=FcorWPL
  )


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

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

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

  #clean up
    rm(sigma, scale)


#-----------------------------------------------------------
#INTEGRAL LENGTH SCALES
  #available variables in REYN$imfl
    vars <- attributes(REYN$imfl)$names

  #scalar length scales
    #variables in imfl
      scal <- c("hor_x", "hor_y", "hor_z", "fst_SONIC_T", "fst_FD_mole_H2O_insitu", "fst_FD_mole_CO2_insitu","FD_mole_CH4_con")
      whr_scal <- sapply(1:length(scal), function(x) which(vars == scal[x]))

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

  #variance length scales
    #calculate scales
      isca_vari <- sapply(whr_scal, function(x) INTsca(
      	  parcel=REYN$data$parcel,
      	  depe=REYN$imfl[,x]^2
      	)
      )
      isca_vari <- as.data.frame( matrix(isca_vari, ncol=length(whr_scal)) )
      attributes(isca_vari)$names <- paste(attributes(REYN$data)$names[whr_scal], "_var", sep="")

  #flux length scales
    #variables in imfl
      whr_flux <- c("Ht", "FD_mole_CH4_flux", "FD_mole_H2O_flux", "FD_mole_CO2_flux", "E")
      whr_flux <- sapply(1:length(whr_flux), function(x) which(vars == whr_flux[x]))
  
    #calculate scales
      isca_flux <- sapply(whr_flux, function(x) INTsca(
          parcel=REYN$data$parcel,
          depe=REYN$imfl[,x]
        )
      )
      isca_flux <- as.data.frame( matrix(isca_flux, ncol=length(whr_flux)) )
      attributes(isca_flux)$names <- attributes(REYN$data)$names[whr_flux]

  #assign to list
    REYN$isca <- as.data.frame(c(isca_scal, isca_vari, isca_flux))

  #clean up
    rm(vars, scal, whr_scal, whr_flux, isca_scal, isca_vari, isca_flux)


#-----------------------------------------------------------
#STATISTICAL ERROR FOR FLUXES
  REYN$erro <- REYNerro(
    isca=REYN$isca,		#integral scale lengths
    mn=REYN$mn,			#mean values
    rc=REYN$cor,		#flux correlation coefficient
    Lf=REYN$max$parcel,		#averaging distance (length of air parcel flown through [m])
    wplc=REYN$wplc,
    FcorWPL=FcorWPL
  )


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


}