# attributes(data)$names # [1] "t_utc" "d_x_utm" "d_y_utm" "d_z_m" "d_xy_travel" "d_xy_flow" "PSI_aircraft" "uvw" # [9] "u_met" "v_met" "w_met" "p_air" "T_air" "FD_mole_H2O" "FD_mole_CH4" "d_z_terrain" # [17] "d_z_ABL" "T_surface" "R_SW_down" REYNflux_FD_mole_dry <- function( data=eddy.data, basetype=c("mean", "trend", "3order")[1], FcorPOT=TRUE, FcorPOTl=NULL ) { #calculate the accumulated length of the parcel of air that has moved past the tower #calculate horizontal wind speed dum_parcel <- sqrt(data$fst_u^2 + data$fst_v^2) #accumulated parcel data$parcel <- cumsum(dum_parcel/freq) ############################################################ #THERMODYNAMICS: MOISTURE, DENSITIES AND TEMPERATURES ############################################################ #----------------------------------------------------------- #GENERAL CONVERSIONS #partial pressure of water vapor data$p_H2O <- data$slow_p * data$fst_FD_mole_H2O_insitu / (1 + data$fst_FD_mole_H2O_insitu) #partial density of H2O data$rho_H2O <- data$p_H2O / R_igl / data$slow_Temp #total (wet) air density [mol m-3] data$rho_air <- data$slow_p / R_igl / data$slow_Temp #dry air density [mol m-3] data$rho_dry <- data$rho_air - data$rho_H2O #virtual temperature -> the temperature of a moist air parcel at which a theoretical #dry air parcel would have a total pressure and density equal to the moist parcel of air [K] #http://en.wikipedia.org/wiki/Virtual_temperature # data$T_v <- data$Temp * (1 + 0.61 * data$q) data$T_v <- data$slow_Temp / (1 - ((data$p_H2O / data$slow_p) * (1 - mvmd)) ) #----------------------------------------------------------- #CONSIDER HUMIDITY IN DRY ADIABATIC CONSTANT #cp as function of humidity (Webb 1980 Eq 40) [J kg-1 K-1] == [m2 s-2 K-1] #dry mole fraction fomulation identical to within 1e-13 #data$cph <- (cpd * data$rho_dry + cpv * data$rho_H2O) / data$rho_air data$cph <- (cpd * 1 + cpv * data$fst_FD_mole_H2O_insitu) / (1 + data$fst_FD_mole_H2O_insitu) #cv as function of humidity [J kg-1 K-1] == [m2 s-2 K-1] #data$cvh <- (cvd * data$rho_dry + cvv * data$rho_H2O) / data$rho_air data$cvh <- (cvd * 1 + cvv * data$fst_FD_mole_H2O_insitu) / (1 + data$fst_FD_mole_H2O_insitu) #specific gas constant for humid air [J kg-1 K-1] == [m2 s-2 K-1] data$Rh <- data$cph - data$cvh #Kappa exponent for ideal gas law (Poisson) [-] data$Kah <- data$Rh / data$cph mn <- data.frame(Kah=mean(data$Kah, na.rm=TRUE)) #----------------------------------------------------------- #POTENTIAL TEMPERATURE AND DENSITIES #potential temperature at 1000 hPa [K] data$T_air_0 <- poisP2T(T1=data$slow_Temp, p1=data$slow_p, p2=p0, kappa=mn$Kah) #virtual potential temperature at 1000 hPa [K] data$T_v_0 <- poisP2T(T1=data$T_v, p1=data$slow_p, p2=p0, kappa=mn$Kah) #use potential temperature and densities? if(FcorPOT == TRUE) { #define pressure level plevel <- ifelse(!is.null(FcorPOTl), FcorPOTl, mean(data$slow_p, na.rm=TRUE) ) #potential temperature at mean pressure level data$slow_Temp <- poisP2T(T1=data$slow_Temp, p1=data$slow_p, p2=plevel, kappa=mn$Kah) #potential densities at mean pressure level #dry air data$rho_dry <- poisP2rho(rho1=data$rho_dry, p1=data$slow_p, p2=plevel, kappa=mn$Kah) #H2O data$rho_H2O <- poisP2rho(rho1=data$rho_H2O, p1=data$slow_p, p2=plevel, kappa=mn$Kah) #wet air data$rho_air <- poisP2rho(rho1=data$rho_air, p1=data$slow_p, p2=plevel, kappa=mn$Kah) #clean up rm(plevel) } ############################################################ #TIME SERIES AVERAGES ############################################################ #data frame mn <- as.data.frame( matrix(colMeans(data, na.rm=TRUE), ncol=ncol(data)) ) attributes(mn)$names <- attributes(data)$names #aircraft heading as vector average mn$PSI_aircraft <- vec2az(matrix(colMeans(az2vec(data$PSI_aircraft*d2r), na.rm=TRUE), ncol=2)) * r2d #wind direction as vector average data$PSI_uv <- vec2az(matrix(c(data$v_met, data$u_met), ncol=2)) mn$PSI_uv <- vec2az(matrix(c(mn$v_met, mn$u_met), ncol=2)) ############################################################ #ROTATION INTO THE MEAN WIND ############################################################ #rotation angle rotang <- ((mn$PSI_uv +180 ) * d2r) %% (2*pi) #transformation matrix B <- matrix(nrow=3, ncol=3) B[1,1] <- cos(rotang) B[1,2] <- sin(rotang) B[1,3] <- 0. B[2,1] <- -sin(rotang) B[2,2] <- cos(rotang) B[2,3] <- 0. B[3,1] <- 0. B[3,2] <- 0. B[3,3] <- 1. BT <- t(B) #wind in (horizontal) geodetic coordinates U <- rbind(data$v_met, data$u_met, data$w_met) #actual rotation Urot <- B %*% U data$u_hor <- Urot[1,] data$v_hor <- Urot[2,] data$w_hor <- Urot[3,] mn$u_hor <- mean(Urot[1,], na.rm=TRUE) mn$v_hor <- mean(Urot[2,], na.rm=TRUE) mn$w_hor <- mean(Urot[3,], na.rm=TRUE) rm(rotang, U, Urot) ############################################################ #BASE STATE AND DEVIATIONS ############################################################ #base state #basetype <- c("mean", "trend", "3order")[2] base <- sapply(1:ncol(data), function(x) base.state(data$parcel, data[,x], basetype)) base <- as.data.frame(matrix(base, ncol=ncol(data))) attributes(base)$names <- attributes(data)$names #vector averages for azimuth angles if basetype == "mean" if(basetype == "mean") { base$PSI_aircraft <- mn$PSI_aircraft base$PSI_uv <- mn$PSI_uv } #str(base) #immidiate fluctuations imfl <- sapply(1:ncol(data), function(x) data[,x] - base[,x]) imfl <- as.data.frame(matrix(imfl, ncol=ncol(data))) attributes(imfl)$names <- attributes(data)$names #correct wind direction from (detrended) wind components PSI_uv_dum <- vec2az(matrix(c(imfl$v_met + mn$v_met, imfl$u_met + mn$u_met), ncol=2)) imfl$PSI_uv <- (PSI_uv_dum - mn$PSI_uv) rm(PSI_uv_dum) #same should be done for PSI_aircraft #variances (corresponding to base state treatment) sd <- as.data.frame( matrix(sqrt(colVars(imfl, na.rm=TRUE)), ncol=ncol(imfl)) ) attributes(sd)$names <- attributes(imfl)$names ############################################################ #ROTATION OF STRESS TENSOR ############################################################ #wind deviations in MET coordinates dx <- imfl$v_met dy <- imfl$u_met dz <- imfl$w_met #stress tensor M <- rbind( c(mean(dx * dx, na.rm = TRUE), mean(dx * dy, na.rm = TRUE), mean(dx * dz, na.rm = TRUE)), c(mean(dy * dx, na.rm = TRUE), mean(dy * dy, na.rm = TRUE), mean(dy * dz, na.rm = TRUE)), c(mean(dz * dx, na.rm = TRUE), mean(dz * dy, na.rm = TRUE), mean(dz * dz, na.rm = TRUE)) ) #rotation into mean wind coordinate system Mrot1 <- B %*% M Mrot2 <- Mrot1 %*% BT #clean up rm(B, BT, dx, dy, dz, M, Mrot1) ############################################################ #MOMENTUM FLUX AND FRICTION VELOCITY ############################################################ #----------------------------------------------------------- #FROM INITIAL COMPONENTS #immidiate fluxes from deviations in mean wind coordinates imfl$u_star2_x <- -(imfl$u_hor * imfl$w_hor) imfl$u_star2_y <- -(imfl$v_hor * imfl$w_hor) imfl$u_star <- NaN #----------------------------------------------------------- #FROM STRESS TENSOR #u_star [m s-1]; optionally only considers the along wind stress u_star_x; Foken (2008) Eq.(2.23) mn$u_star2_x <- -Mrot2[1,3] mn$u_star2_y <- -Mrot2[2,3] mn$u_star <- (mn$u_star2_x^2 + mn$u_star2_y^2)^(1/4) #TK3 style, necessary if flight path not with wind #correlations cor <- data.frame( u_star2_x= -Mrot2[1,3] / sd$u_hor / sd$w_hor, u_star2_y= -Mrot2[2,3] / sd$v_hor / sd$w_hor ) #wind variance; deviations from initials sd are < 2% sd_dum <- sqrt(abs(diag(Mrot2))) sd$u_hor <- sd_dum[1] sd$v_hor <- sd_dum[2] sd$w_hor <- sd_dum[3] #----------------------------------------------------------- #CLEAN UP rm(Mrot2, sd_dum) ############################################################ #SENSIBLE HEAT FLUX ############################################################ #SENSIBLE HEAT FLUX #flux in kinematic units [K m s-1] imfl$F_H_kin <- imfl$w_hor * imfl$slow_Temp mn$F_H_kin <- mean(imfl$F_H_kin, na.rm=TRUE) #conversion to units of energy [W m-2] == [kg s-3] imfl$F_H_en <- (cpd * base$rho_dry * Md + cpv * base$rho_H2O * Mv) * 1e-3 * imfl$F_H_kin mn$F_H_en <- mean(imfl$F_H_en, na.rm=TRUE) #BUOYANCY FLUX considering water vapor buoyancy and 1000 hPa reference (virt. pot. temp.) -> z/L #flux in kinematic units [K m s-1] imfl$F_H_kin_v_0 <- imfl$w_hor * imfl$T_v_0 mn$F_H_kin_v_0 <- mean(imfl$F_H_kin_v_0, na.rm=TRUE) #CORRELATIONS cor$F_H_kin <- stats::cor(imfl$w_hor, imfl$slow_Temp, use="pairwise.complete.obs") cor$F_H_en <- cor$F_H_kin cor$F_H_kin_v_0 <- stats::cor(imfl$w_hor, imfl$T_v_0, use="pairwise.complete.obs") ############################################################ #LATENT HEAT FLUX ############################################################ #latent heat flux in kinematic units [mol m-2 s-1] imfl$F_LE_kin <- base$rho_dry * imfl$w_hor * imfl$fst_FD_mole_H2O_insitu mn$F_LE_kin <- mean(imfl$F_LE_kin, na.rm=TRUE) #latent heat of vaporization (Eq 2.55 Foken 2008) [J kg-1] == [m2 s-2] Lv <- 2500827 - 2360 * (base$slow_Temp - T0) #latent heat flux in units of energy [W m-2] == [kg s-3] imfl$F_LE_en <- Lv * Mv * 1e-3 * imfl$F_LE_kin mn$F_LE_en <- mean(imfl$F_LE_en, na.rm=TRUE) #correlation cor$F_LE_kin <- stats::cor(imfl$w_hor, imfl$fst_FD_mole_H2O_insitu, use="pairwise.complete.obs") cor$F_LE_en <- cor$F_LE_kin #clean up rm(Lv) ############################################################ #CH4 FLUX ############################################################ #CH4 flux in kinematic units [mol m-2 s-1] imfl$F_CH4_kin <- base$rho_dry * imfl$w_hor * imfl$FD_mole_CH4 mn$F_CH4_kin <- mean(imfl$F_CH4_kin, na.rm=TRUE) #CH4 flux in mass units [mg m-2 h-1] imfl$F_CH4_mass <- imfl$F_CH4_kin * M_CH4 * 1e3 * 3600 mn$F_CH4_mass <- mean(imfl$F_CH4_mass, na.rm=TRUE) #correlation cor$F_CH4_kin <- stats::cor(imfl$w_hor, imfl$FD_mole_CH4, use="pairwise.complete.obs") cor$F_CH4_mass <- cor$F_CH4_kin ############################################################ #AUXILARY PARAMETERS ############################################################ #turbulence intensity from total wind vector (Stull, Eq. 1.4d) #total wind vector tot <- sqrt(data$u_met^2 + data$v_met^2 + data$w_met^2) #turbulence intensity should be <0.5 to allow for Taylors hypothesis mn$I <- sd(tot, na.rm=TRUE) / mean(tot, na.rm=TRUE) #Obukhov length (used positive g!) mn$d_L_v_0 <- (-( ( (mn$u_star)^3 / ( kap * g / mn$T_v_0 * mn$F_H_kin_v_0 ) ) )) #stability mn$sigma <- mn$uls_z / mn$d_L_v_0 #convective (Deardorff) velocity [m s-1] #missing values in Deardorff velocity and resulting variables when buoyancy flux is negative! mn$w_star <- ( g * mn$d_z_ABL / mn$T_v_0 * mn$F_H_kin_v_0 )^(1/3) #(free) convective time scale [s], often in the order of 5-15 min mn$t_star <- mn$d_z_ABL / mn$w_star #temperature scale (eddy temperature fluctuations) [K] #surface layer #mn$T_star_SL <- - mn$F_H_kin_v_0 / mn$u_star #according to Stull (1988) p. 356 mn$T_star_SL <- - mn$F_H_kin / mn$u_star #according to Foken (2008) p.42, fits with ITC assessment #mixed layer mn$T_star_ML <- mn$F_H_kin / mn$w_star #according to Stull (1988) p. 356 #humidity scale (eddy moisture fluctuations) [mol mol-1 dry air] #surface layer mn$FD_mole_H2O_star_SL <- - mean(mn$F_LE_kin / base$rho_dry, na.rm=TRUE) / mn$u_star #mixed layer layer mn$FD_mole_H2O_star_ML <- mean(mn$F_LE_kin / base$rho_dry, na.rm=TRUE) / mn$w_star #clean up rm(base, tot) ############################################################ #EXPORT RESULTS ############################################################ #PREPARE DATA #mins mi <- as.data.frame( matrix(colMins(data, na.rm=TRUE), ncol=ncol(data)) ) attributes(mi)$names <- attributes(data)$names #maxs ma <- as.data.frame( matrix(colMaxs(data, na.rm=TRUE), ncol=ncol(data)) ) attributes(ma)$names <- attributes(data)$names #assemble export list export <- list( data=data, #data including internal calculations min=mi, #min max=ma, #max mn=mn, #mean sd=sd, #standard deviation imfl=imfl, #immidiate fluctuations cor=cor #correlation coefficient ) #clean up rm(data, mi, ma, mn, sd, imfl, cor) #return result return(export) }