extract <- function(x, y, data, xmin, xmax, ymin, ymax, t=1) {

# Extract a subset of data from either a two-dimensional map (matrix) with p
# longitude and q latitude points, respectively, or a three-dimensional
# array (field) with p longitude and q latitude as the first two space 
# dimensions and n (time) as the third dimension, by providing the latitude and
# longitude range of interest and a subset of time slices 
#
# Description:
#
#      Returns a subset of a two-dimensional map (matrix) or a subset
#      of a three-dimensinal array (field) 
#
# Usage:
#
#      extract(x, y, data, xmin, xmax, ymin, ymax, t=1)
#
# Arguments:
#
#      x: Vector of longitudes 
#
#      y: Vector of latitudes
#
#      data: Two-dimensional map (matrix) or three-dimensional array of data.
#            Space dimensions of `data' must be `length(x)' by `length(y)'
#            respectively. 
#            Three-dimensional array must have the following dimension order:
#            first dimension for longitude, second dimension for latitude,
#            and third dimension for time
#
#      xmin: Lower bound longitude from where to start extracting data
#
#      xmax: Upper bound longitude up to where to extract data
#
#      ymin: Lower bound latitude from where to start extracting data
#
#      ymax: Upper bound latitude up to where to extract data
#
#      t: Subset of times slices to be extracted. Default is 1 (i.e. extract
#         only the first time slice
#
# Output:
#
#      A list with four components is returned invisibly:
#
#         x: Vector of longitudes within the selected range
#
#         y: Vector of latitudes within the selected range
#
#         data: Two-dimensional (matrix) or three-dimensional (array)
#               containing the selected data
#
#         t: Extracted times slices 
#
#
# Author:
#       
#      Caio Coelho <c.a.d.s.coelho@reading.ac.uk> 28 October 2005
#      Chris Ferro <c.a.t.ferro@reading.ac.uk> 
#
# Examples:
#
#      # Extract data from a three-dimensional array (field)
#      x <- seq(-20, 20, 5)
#      y <- seq(30, 60, 5)
#      dim <- c(length(x), length(y), 100)
#      data <- array(rnorm(prod(dim)), dim)
#      output <- extract(x, y,data,xmin=-12,xmax=12,ymin=32,ymax=58,t=c(1:3))
#      output <- extract(x, y,data,xmin=-12,xmax=6,ymin=32,ymax=50,t=c(2,4,6))
#      output <- extract(x, y,data,xmin=-18,xmax=11,ymin=36,ymax=59,t=c(10,20,30,40))
#
#      # Extract data from a two-dimensional map (matrix)
#      matr <- data[,,1]
#      output <- extract(x, y,matr,xmin=-12,xmax=12,ymin=32,ymax=58)
#      output <- extract(x, y,matr,xmin=-12,xmax=6,ymin=32,ymax=50)


if (length(dim(data))==2) {
  xnew <- x
  ynew <- y 
    xnew <- x[x >= xmin & x <= xmax]
    ynew <- y[y >= ymin & y <= ymax]
    datanew <- data[match(xnew, x), match(ynew, y), drop = FALSE] 
}  

if(length(dim(data))==3) {
  xnew <- x
  ynew <- y
  tnew <- t
    xnew <- x[x >= xmin & x <= xmax]
    ynew <- y[y >= ymin & y <= ymax]
    datanew <- data[match(xnew, x), match(ynew, y), t, drop = FALSE]
} 

invisible(list(x = xnew, y = ynew, data = datanew, t = t))
}