Silicate, phosphate and nitrate

## ====================================================================================
## Nutrient concentrations at surface
## The world Ocean Database
## implemented by Karline Soetaert
## ====================================================================================
library(OceanView)

## Load if not yet loaded
library(RNetCDF)
Directory <- "ftp://ftp.nodc.noaa.gov/pub/data.nodc/woa/WOA05nc/annual/"
url <- Directory
ftplist <- c("i00an1", "n00an1", "O00an1", 
             "p00an1", "s00an1", "t00an1")

path <- getwd()
Listfiles <- list.files(path)
for (i in 1:length(ftplist)) {
  File  <- paste(ftplist[i], ".nc", sep = "") 
  if (!File %in% Listfiles) {
    Url       <- paste(Directory, File, sep = "")
    download.file(Url, File, method =  "auto", quiet = FALSE, 
      mode = "wb", cacheOK = TRUE)   
  }
}

# read data and put in list

readWOA <- function(input) {
  D.nc <- open.nc(paste(path,"/",input, ".nc", sep = ""))
  lon      <- var.get.nc(D.nc, 'lon')
  lat      <- var.get.nc(D.nc, 'lat')     
  time     <- var.get.nc(D.nc, 'time')     
  depth    <- var.get.nc(D.nc, 'depth')     
  value    <- var.get.nc(D.nc, input)
  name     <- att.get.nc(D.nc, input, "long_name")
  units    <- att.get.nc(D.nc, input, "units")
 # create list with data
  return (list(lon = lon, lat = lat, depth = depth, time= time,
      value = value, name = name, units = units))
}
silicate  <- readWOA("i00an1") 
nitrate   <- readWOA("n00an1") 
oxygen    <- readWOA("O00an1") 
phosphate <- readWOA("p00an1")

# plot

library(plot3D)
par(mfrow = c(2, 2), oma = c(2, 0, 0, 0))
with (silicate, image2D(z = value[ , ,1], x = lon, y = lat, main = name, 
      xlab = "longitude", ylab = "latitude", clab = c("", "", units)))
with (phosphate, image2D(z = value[ , ,1], x = lon, y = lat, main = name, 
      xlab = "longitude", ylab = "latitude", clab = c("", "", units)))
with (nitrate, image2D(z = value[ , ,1], x = lon, y = lat, main = name, 
      xlab = "longitude", ylab = "latitude", clab = c("", "", units)))
abline(v = 180, lwd = 2)
with (nitrate, image2D(z = value[180 , , ], x = lat, y = depth, 
        ylim = c(5500, 0), NAcol = "black", resfac = 2,
        xlab = "dgN", ylab = "depth, m",
        main = c("nitrate", paste(lon[180], "dgE")), 
        clab = c("", "", units)))