library(ncdf4) 		## for NetCDF files
library(RColorBrewer)   ## for some nice color scales, see http://colorbrewer2.org
library(plyr)		## for the aaply family of functions

source("r/image.plot.r") ## for lat/lon plots

## a useful polynomial approximation for saturation water vapor pressure as a function of T
e.sat.water <- function(T, p) {
  T.C <- T - 273.15
  e.sat <- 6.107799961 + T.C*(4.436518521e-1 + T.C*(1.428945805e-2 + T.C*(2.650648471e-4 + T.C*(3.031240396e-6 + T.C*(2.034080948e-8 + 6.136820929e-11*T.C))))) ## polynomial approximation taken from Appendix 4, Pruppacher & Klett, 2nd ed; valid for -50 deg C < T < 50 deg C
}

## ERA-Interim reanalysis temperature field
nc.era <- nc_open("~quaas/data/ERA__Interim__T_GDS0_ISBL_123__1.5x1.5xL37__198901-200712.nc")
lev.era <- ncvar_get(nc.era, "lv_ISBL1")
lon.era <- ncvar_get(nc.era, "g0_lon_3")
lat.era <- ncvar_get(nc.era, "g0_lat_2")
time.era <- ncvar_get(nc.era, "time") * 3600 + as.POSIXlt("1800-01-01", tz = "UTC") ## translate from the strange ERA units
T.era <- ncvar_get(nc.era, "T_GDS0_ISBL_123")

print(dim(T.era))  ## 240 lon x 121 lat x 37 pressure levels x 228 months... a big field

## plot a vertical temperature and potential temperature profile in the tropics (120 deg W on the equator)
plot(T.era[which(lon.era == -120), which(lat.era == 0), , 1], lev.era, type = "l", col = "red",
     ylim = c(1000, 0), xlim = c(200, 350), xlab = "K", ylab = "p (hPa)")
lines(T.era[which(lon.era == -120), which(lat.era == 0), , 1] * (1000 / lev.era)^(2 / 7), lev.era, col = "blue")
legend("bottomleft", title = "Tropical ocean", c("temperature", "potential temperature"), col = c("red", "blue"), lty = "solid")

## look at the 1000 hPa temperature (near-surface temperature) in June, July and August (JJA) and December, January, February (DJF)
Tsurf.era.JJA <- T.era[,, which(lev.era == 1000), months(time.era) %in% c("June", "July", "August")]
Tsurf.era.DJF <- T.era[,, which(lev.era == 1000), months(time.era) %in% c("December", "January", "February")]

## calculate the temporal mean of the JJA and DJF near-surface temperature fields
Tsurf.era.JJA.bar <- apply(Tsurf.era.JJA, c(1,2), mean)
Tsurf.era.DJF.bar <- apply(Tsurf.era.DJF, c(1,2), mean)

image.map(lon.era, lat.era, Tsurf.era.JJA.bar, col = brewer.pal(5, "OrRd"), zlim = c(220, 320))
image.map(lon.era, lat.era, Tsurf.era.DJF.bar, col = brewer.pal(5, "OrRd"), zlim = c(220, 320))

## calculate the temporal mean of the temperature field
T.era.bar <- apply(T.era, c(1,2,3), mean)
## calculate the zonal mean of the temporal-mean temperature field
T.era.bar.bracket <- apply(T.era.bar, c(2,3), mean)

## calculate zonal-temporal-mean potential temperature 
theta <- aaply(T.era.bar.bracket, 1, function(T) T * (1000 / lev.era)^(2 / 7))

## plot
contour(lat.era, lev.era, T.era.bar.bracket, ylim = c(1000, 0), xlab = "Latitude", ylab = "p (hPa)")

## plot, this time into a PDF file:
pdf("T-contour.pdf", width = 8, height = 5)
contour(lat.era, lev.era, T.era.bar.bracket, ylim = c(1000, 0), xlab = "Latitude", ylab = "p (hPa)")
dev.off() ## close the output file