[This document is not complete -- it only documents the turbulent flux dat routines at present. If users find it useful, more can be added.]

Introduction

NCAR/RTF/Surface Group processes most of our data using a commercial data analysis package, S+, or its open-source look-alike, R. We have built an extensive library of routines to work with our data using these packages, which are now based on a object-oriented model. We have defined a "dat" object, which has several attributes that are available for our routines to query:

We also have defined a more primative "ats" object without the station attribute or "dat" class. "dat" objects usually are either read from a NetCDF file or calculated from other "dat" objects. We have developed functions to allow a set of "dat" objects that are derived to be referenced as though they had been read in. For example, dat("T") reads air temperature directly from the NetCDF file, but dat("Q") computes specific humidity from dat("T"), dat("RH"), and dat("P"), where dat("RH") is the relative humidity measured by our sensors and dat("P") is the pressure also measured. dat("Q") is created by calling a routine "dat.Q" which we have written.

A description of some of our "dat.X" functions is given in this document. Throughout this document, code which is only "housekeeping" in nature has been deleted to focus on the calculations themselves.

Parameters

Calculation by the "dat" functions is controlled by a set of parameters set by the "dpar()" function or a point-and-click equivalent "dpopup()". Through these parameters, the time period of interest (usually beginning time and duration) is specified. Additional options can select the measurement location (station and/or height), averaging period (see next section).

Of particular interest to the below discussion is the "robust" parameter. Many of the routines below require several data inputs to implement a full correction to the data (the dat.Q example above requires 3 input measurements). However, if any of these measurements are missing, the output of a "dat" routine necessarily must be NA. For some applications (such as when we are in the field), it is useful to get a "quick-and-dirty" answer whether or not all the measurements are available. "robust=T" tells the dat functions to compute a result based on the fewest inputs. "robust=F" forces the dat functions to perform all data corrections, but may yield fewer results.

Averaging

If "avg" is specified by dpar() or dpopup(), the dat routines will compute a true average of the values from values stored in the NetCDF files. All of our NetCDF files to date were created using 5 minute averages. Thus dat("T") with avg=1800 [seconds] would read data from 6 5-minute blocks and create a 30-minute average.

The averaging process is smart enough to include the "weights" parameter so that if only 100 samples were averaged in one 5 minute block and 300 samples in another, the first would be weighted only 1/3 as much.

If the dpar parameter "simple.avg" or "robust" is "F", the average of higher- order momements (variance, covariance, skewness, etc.) will include the correct variation of the 5-minute average values in the computation.

Air Density

rho.air

Many routines require the density of air, and we usually measure temperature, humidity, and pressure, so that this can be computed. If robust=T, only one (median) value is returned to keep matching simple.

dat.rho.air <- function(x,cache=F,robust=dpar("robust"),...)
{
  # density of (moist) air (kg/m^3)
  R <- 287		# J/kg-K
  p <- dat("P")*100	# pascals
  q <- dat("Q")
  Td <- dat("T")
  Tv <- (273.15 + Td)*(1 + 0.608*q*1e-3)
  x <- p/R/Tv
  if (!is.null(robust) && robust) {
    x <- median(x, na.rm=T)
    attr(x,"dunits") <- "kg/m^3"
    }
  else {
    attr(x,"dunits") <- rep("kg/m^3",ncol(x))
  }
  x
}

rho.dry

Sometimes, we only want the density of dry air. As with rho.air, only one value is returned for robust=T.

dat.rho.dry <- function(x,cache=F,robust=dpar("robust"),...)
{
  # density of dry air (kg/m^3)
  R <- 287		# J/kg-K
  mr <- dat("MR")
  p <- dat("P")*100/(1 + dat("MR")/622)		# P - Pv, pascals
  Td <- dat("T")
  TK <- 273.15 + Td
  x <- p/R/TK		
  if (!is.null(robust) && robust) {
    x <- median(x, na.rm=T)
    attr(x,"dunits") <- "kg/m^3"
    }
  else {
    attr(x,"dunits") <- rep("kg/m^3",ncol(x))
  }
  x
}

Turbulent Fluxes

These start getting complicated very fast due to all the corrections and interdependencies. It probably is necessary to read through EVERYTHING in this section to get the complete picture.

u*

This is wrong of course. However, to keep things simple (namely that we always have a u* value and aren't getting imaginary values), we just compute a flux magnitude rather than the component of the stress vector in the streamwise direction even when robust=F. 

"dat.u*" <- function(x,cache=F,...)
{
  x <- (dat("u'w'")^2 + dat("v'w'")^2) ^ (0.25)
  attr(x,"dunits") <- rep("m/s",ncol(x))
  x
}

H

This routine computes the total sensible heat flux, including the heat capacity of moisture in the air parcel. Most of the work in this routine is done by the functions dat.w'h2o', dat.w't', and calc.H, which figure out how to convert the flux of sonic virtual temperature to real temperature, etc.

Missing temperatures (needed to compute the heat capacity) are filled in with the median temperature over the entire time period in the case that robust=T. 

dat.H <- function(x,robust=dpar("robust"),...)
{
  # Turbulent sensible heat flux
  Td <- dat("T")
  wq <- dat("w'h2o'")*1e-3
  wt <- dat("w't'")
  x <- calc.H(wt=wt,wq=wq,Td=Td,rho=dat("rho.air"))
  x
}

dat.H.dry

Some people want to know what the sensible heat would be without the water vapor contribution to the specific heat. As before, missing data are filled in if robust=T. 

dat.H.dry <- function(x,robust=dpar("robust"),...)
{
  # Turbulent sensible heat flux without contribution from varying Cp 
  Q <- dat("Q")
  x <- calc.H(wt=dat("w't'"),Q=dat("Q"),rho=dat("rho.air"))
  x
}

dat.LE

Given all the nonsense below, this is quite simple, with the real work hidden in dat("w'mr'"). Note that LE normally is defined from the specific humidity flux w'q' times dat("rho"), but the approach we use allows us to deal only with one dat routine. 

dat.LE <- function(x,robust=dpar("robust"),...)
{
  # Turbulent latent heat flux
  x <- dat("Lv")*dat("rho.dry")*dat("w'mr'")*1e-3
  attr(x,"dunits") <- rep("W/m^2",ncol(x))
  x
}

w't'

dat.w't' corrects the measured sonic virtual temperature flux for the contribution due to the water vapor flux (Schotanus et al., BLM, 1983). However, since our water vapor flux usually is measured by a krypton hygrometer, it is necessary to apply the Webb correction (Webb et al., QJRMS, 1980) and oxygen correction (Campbell App.C) to the w'h2o' measurements. Since these corrections both need w't' itself, these corrections must be solved for simultaneously. However, the correction can be written out explicitly as a function of w'tc' and w'h2o'. 

"dat.w't'" <- function(x,robust=dpar("robust"),...)
{
  # Corrections to  assume E measured as  (m/s g/m^3)
  if (!is.null(robust) && robust)
    x <- dat("w'tc'")
  else {
    rho.air <- dat("rho.air")
    Q <- dat("Q")*1e-3
    Td <- dat("T")
    TK <- 273.15 + Td
    E <- dat("w'h2o'")*1e-3     # m/s kg/m^3
    o2corr <- calc.o2corr(dat("P"),TK,rho.air,E)
    factor <- 0.51*(1 + Q*(1/0.622-1))
    x <- (dat("w'tc'") - factor*TK*E/rho.air)/(1 + factor*(Q+o2corr))
    x[is.inf(x)] <- NA
    }
  attr(x,"dunits") <- rep("m/s degK",ncol(x))
  x
}

w'mr'

This is the workhorse for dat.LE which solves the second set of correction equations for w'mr' in terms of w'tc' and w'h2o', as done for w't'. 

"dat.w'mr'" <- function(x,cache=F,robust=dpar("robust"),...)
{
  # Assume E measured as  (m/s g/m^3)
  E <- dat("w'h2o'")
  rhod <- dat("rho.dry")
  x <- E/rhod*1e-3
  if (!(!is.null(robust) && robust)) {
    Q <- dat("Q")*1e-3
    MR <- Q/(1-Q)
    Td <- dat("T")
    TK <- 273.15 + Td
    o2corr <- calc.o2corr(dat("P"),TK,dat("rho.air"),E)
    factor <- 0.51*(1 + Q*(1/0.622-1))
    x <- (1 + MR/0.622)*(x + (MR+o2corr/rhod)/TK*dat("w'tc'"))/
	(1 + factor*(Q+o2corr))
    }
  x <- x*1e3		# m/s g/kg
  attr(x,"dunits") <- rep("m/s g/kg",ncol(x))
  x
}

calc.o2corr

This makes the oxygen correction to kryptons (from Campbell, Appendix C and my krypton write-up). I have been told recently that Campbell wrote a later document which reduces the magnitude of this correction.

I have removed the messy logic from this write-up which tests whether kryptons were used (dpar("fh2os")) and skips this routine if not. 

calc.o2corr <- function(P,TK,rho.air,E,robust=dpar("robust")) {
  hs <- dpar("fh2os")
  P <- 100*P
  corr <- (4.7e-5)*P/(rho.air*TK)
  corr
}

calc.H

This deals with the water contribution to heat capacity in the dat.H routines as per Stull eq. 10.7.1b or Stull eq. 10.7.1c and Fleagle and Businger (second edition) eq. 6.23. 

calc.H <- function(wt,wq,Td,Q,rho.air)
{
  # Turbulent sensible heat flux
  # Constant values from Smithsonian Meteorological Tables:
  Cpd <- 1006	# for dry air, range is 1004.8-1007.4 from -50 to +50 C and
                # 700 to 1100 mb, J/kg-K
  Cpv <- 1857	# for water vapor, value for 0 C, range is 1850 to 1870 
                # from -50 to +50C, J/kg-K
  # No fast water covariance
  if (missing(wq)) {
    Hname <- "H.dry"
    Q <- Q/1000   # kg/kg
    # Stull, eq. 10.7.1b
    Cp <- Cpd*(1 + 0.84*Q)
    x <- Cp * rho.air
    x <- x * wt
  }
  else {
    Hname <- "H"
    x <- rho.air
    # Stull, eq. 10.7.1c or Fleagle&Businger (2nd ed) eq. 6.23, 
    # using Cpd&Cpv directly
    x <- x * (Cpd*wt + Cpv*Td*wq)
  }
  dimnames(x) <- list(NULL,paste(Hname,sfxnames(wt,2),sep=""))
  attr(x,"dunits") <- rep("W/m^2",ncol(x))
  x
}

dat.Lv

I note here that the Latent heat of vaporatizaton of water is a weak function of temperature (varies by 1% in 10 C), but appears to be well fit by a linear fit to data I found on the WWW.

However, to avoid matching problems, I'm only using the first temperature data point, which assumes that all the real temperatures aren't TOO much different. 

dat.Lv <- function(x,cache=F,...)
{
  # Still, just 
  if (!is.null(dpar("robust")) & !dpar("robust")) {
  	x <- 2501000 - 2373*dat("T")[,1]
	x[is.na(x)] <- 2501000 - 2373*median(dat("T"),na.rm=T)
	}
  else {
	x <- 2.5e6
	}
  dimnames(x) <- list(NULL,paste("Lv",sfxnames(x,2),sep=""))
  class(x) <- c("Lv",class(x))
  attr(x,"dunits") <- rep("J/kg",ncol(x))
  x
}

References

Campbell, Appendix C

Fleagle, R. and J.A.Businger, An Introduction to Atmospheric Physics, Second Edition, Academic Press

Schotanus, P., 1983, ... Boundary Layer Meteor.

Smithsonian Meteorological Tables

Stull, R. 1988, An Introduction to Boundary Layer Meteorology, Kluwer Academic, Dordrecht, 666 pp.

Webb, E., 1980, Quart. J. Royal Meteor. Soc.,