Hycom temperature

## =============================================================================
## Plotting output of a hydrodynamic model - assume files have been downloaded
## Surface values and 2 transects
## implemented by Karline Soetaert
## =============================================================================

require(OceanView)
require(RNetCDF)

# Read the data - manually downloaded
url <- "ftp://ftp.hycom.org/datasets/ATLd0.08/expt_09.1/data/2d/archv.2006123100_2d.nc" 
Nc <- open.nc("archv.2006123100_2d.nc")
print.nc(Nc)

lat <- var.get.nc(Nc, "Latitude")
lon <- var.get.nc(Nc, "Longitude")
tmix <- var.get.nc(Nc, "mixed_layer_temperature")

Nc3 <- open.nc("archv.2006123100_3zt.nc")
print.nc(Nc3)

depth <- var.get.nc(Nc3, "Depth")
temp  <- var.get.nc(Nc3, "temperature") 

layout(matrix(nrow = 4, ncol = 4, data = c(rep(c(1, 1, 1, 2), 3), 3, 3, 3, 4)))          
pm <- par("mar")

clim <- c(-5, 30)
image2D(x= lon, y = lat, z = tmix, col = jet2.col(100), shade = 0.1,  #rasterImage = TRUE, 
 main = "HYCOM ATLd0.08 temperature", colkey = FALSE, clim = clim,
 xlab="", ylab ="", clab = "m/s", NAcol = grey(0.4))
ix <- 902
iy <- 902
abline(h = lat[iy])
abline(v = lon[ix])

# transect 
par(mar = pm + c(0, 0, -2, 0))
image2D(temp[,iy,], col = jet2.col(100), x = lon, y = depth, 
  ylim = rev(range(depth)), xlab = "", ylab  = "",
  main = expression(paste("39.4" , degree*N)), 
    NAcol = grey(0.4), resfac = 2, clim = clim, colkey = FALSE)

par(mar = pm + c(0, -2, 0, 0))
image2D(t(temp[ix,,]), col = jet2.col(100), y = lat, x = depth, 
  xlim = rev(range(depth)),  xlab = "", ylab  = "", 
  main = expression(paste("25.9" , degree*W)), , 
    NAcol = grey(0.4), resfac = 2, clim = clim, colkey = FALSE)

par(mar = pm)
colkey(clab = expression(degree*C), clim = clim, length = 0.8)