R/source.roms.r

In tbrycekelly/TheSource: TheSource: It Will Hold Your Beer

## Set of useful R functions for working with Regional Ocean
## Model System (ROMS) data and netcdf files.
##
## Including LTRANS files.
##
## Author: Thomas Bryce Kelly (tbk14 at fsu.edu)
## http://about.tkelly.org/
##
## Dept of Earth, Ocean & Atmospherical Sciences
## Florida State University
##
## Center for Ocean & Atmospheric Prediction Studies
## Florida State University
##
## National High Magnetic Field Laboratory
## Florida State University



##########################
##### LOADING FUNCTIONS
##########################

#' @title Load in ROMS physics
#' @export
#' @author Thomas Bryce Kelly
#' @import ncdf4
load.roms.physics = function(path, verbose = T) {

    model = load.nc(file = path,
                    var = c('u', 'Huon', 'v', 'Hvom', 'w', 'AKt', 'Akt_bak',
                            'temp', 'rho', 'salt', 'hice', 'f', 'nl_tnu2', 'snow_thick',
                            'uice', 'vice', 'zeta', 'ocean_time'),
                    verbose = verbose)

    model$meta = list(
        path = path,
        time = Sys.time(),
        Source.version = packageVersion('TheSource'),
        R.version = R.version.string
    )
    ## Return
    model
}

#' @title Load in ROMS grid
#' @export
#' @author Thomas Bryce Kelly
#' @import ncdf4
load.roms.grid = function(path, verbose = T) {

    model = load.nc(file = path,
                    var = c('h', 'hc', 'theta_s', 'theta_b',
                            'lat_psi', 'x_psi', 'lon_psi', 'y_psi',
                            'zice', 'Cs_r', 'Cs_w',
                            'mask_psi', 'mask_rho', 'mask_u', 'mask_v',
                            'pm', 'pn', 's_rho', 's_w',
                            'x_psi', 'x_rho', 'x_u', 'x_v',
                            'y_psi', 'y_rho', 'y_u', 'y_v',
                            'Vtransform', 'Vstretching'),
                    verbose =  verbose)

    model$meta = list(
        path = path,
        time = Sys.time(),
        Source.version = packageVersion('TheSource'),
        R.version = R.version.string
    )
    ## Return
    model
}



#' @title Load in ROMS temporal info
#' @export
#' @author Thomas Bryce Kelly
#' @import ncdf4
load.roms.time = function(path, verbose = T) {

    model = load.nc(file = path,
                    var = c('dstart', 'dt', 'dtfast',
                            'nAVG', 'ndefAVG', 'ndefHIS', 'ndtfast',
                            'nHIS', 'nRST', 'ntimes', 'ntsAVG', 'ocean_time'),
                    verbose =  verbose)

    model$meta = list(
        path = path,
        time = Sys.time(),
        Source.version = packageVersion('TheSource'),
        R.version = R.version.string
    )
    ## Return
    model
}



###########################
#### CONVERSION FUNCTIONS
###########################

#' @title Convert ROMS Time
#' @author Thomas Bryce Kelly
#' @param x The datetime numeric to be converted
#' @param tz The timezone for the conversion, almost always UTC
#' @export
conv.time.roms = function(x, origin = "1900-01-01", tz = 'UTC') {
    as.POSIXct(x, origin = origin, tz = tz)
}


#' @title Get ROMS Time
#' @export
#' @author Thomas Bryce Kelly
calc.roms.time = function(time.data, tz = 'UTC') {

    ## If we have ocean_time, then use it, if not try and calculate it based on time steps.
    if (!is.na(time.data$ocean_time[1])) {
        time = conv.time.roms(time.data$ocean_time)
    } else {
        nhist = time.data$nHIS # dt * nHist = time between records
        dstart = time.data$dstart # days since 1900-01-01 00:00:00
        dt = time.data$dt # Sec
        ntimes = time.data$ntimes # number of time frames


        ## Calculate the times (seconds from 1900-01-01)
        time = seq(from = dstart, by = dt * nhist / 86400, length.out = ntimes/nhist) * 86400
    }

    time
}


#' @export
calc.roms.depths = function(grid.data, zeta = NULL, verbose = T) {

    ## Setup
    nx = dim(grid.data$h)[1]
    ny = dim(grid.data$h)[2]
    h = destagger.grid(grid.data$h)
    if (!is.na(grid.data$zice[1])) {
        zice = destagger.grid(grid.data$zice)
    } else {
        zice = matrix(0, nrow = nx - 1, ncol = ny - 1)
    }


    ## Free surface
    if (is.null(zeta)) {
        zeta = matrix(0, nrow = nx - 1, ncol = ny - 1)
    } else {
        if (length(dim(zeta)) > 2) {
            message(' Zeta has dimensions greater than 2. Ignoring zeta...')
            zeta = matrix(0, nrow = nx - 1, ncol = ny - 1)
        } else {
            zeta = destagger.grid(zeta)
        }
    }

    ## Setup depths and calculate
    z = array(0, dim = c(nx-1, ny-1, length(grid.data$Cs_r)))

    if (is.na(grid.data$Vtransform) | grid.data$Vtransform == 1) { ## Equation 1
        if (is.na(grid.data$Vtransform)) {message(' No Vtransform equation found, assuming equation 1!')}
        for (i in 1:length(grid.data$Cs_r)) {
            z0 = (grid.data$s_rho[i] - grid.data$Cs_r[i]) * grid.data$hc + grid.data$Cs_r[i] * (h + zice)
            z[,,i] = z0 + zeta * (1 + z0 / (h + zice)) + zice
        }
    } else if (!is.na(grid.data$Vtransform) & grid.data$Vtransform == 2) { ## Equation 2
        for (i in 1:length(grid.data$Cs_r)) {
            z0 = (grid.data$hc * grid.data$s_rho[i] + grid.data$Cs_r[i]* (h + zice)) / (grid.data$hc + h + zice)
            z[,,i] = zeta + (zeta + h + zice) * z0
        }
    } else {
        message(' No valid Vtransform given. Returning zeros!')
    }

    z
}


#' @export
calc.roms.areas = function(grid.data, zeta = NULL, verbose = T) {

    ## Setup
    nx = dim(grid.data$h)[1]
    ny = dim(grid.data$h)[2]

    if (verbose) { message('Calculating areas for all planes.')}
    h.center = destagger.grid(grid.data$h)
    h.x = destagger.grid.x(grid.data$h)
    h.y = destagger.grid.y(grid.data$h)

    if (!is.na(grid.data$zice[1])) {
        zice = destagger.grid(grid.data$zice)
        zice.x = destagger.grid.x(grid.data$zice)
        zice.y = destagger.grid.y(grid.data$zice)
    } else {
        zice = matrix(0, nrow = nx - 1, ncol = ny - 1)
        zice.x = matrix(0, nrow = nx - 1, ncol = ny)
        zice.y = matrix(0, nrow = nx, ncol = ny - 1)
    }


    ## Free surface
    if (is.null(zeta)) {
        zeta = matrix(0, nrow = nx - 1, ncol = ny - 1)
        zeta.x = matrix(0, nrow = nx - 1, ncol = ny)
        zeta.y = matrix(0, nrow = nx, ncol = ny - 1)

    } else {
        if (verbose) { message(' Calculating free surface...')}
        if (length(dim(zeta)) > 2) {
            message(' Zeta has dimensions greater than 2. Ignoring zeta...')
            zeta = matrix(0, nrow = nx - 1, ncol = ny - 1)
        } else {
            zeta.x = destagger.grid.x(zeta)
            zeta.y = destagger.grid.y(zeta)
            zeta = destagger.grid(zeta)
        }
    }

    ## Setup depths and calculate
    z = array(0, dim = c(nx-1, ny-1, length(grid.data$Cs_r)))
    z.x = array(0, dim = c(nx-1, ny, length(grid.data$Cs_r)))
    z.y = array(0, dim = c(nx, ny-1, length(grid.data$Cs_r)))

    area.f = array(0, dim = c(nx-1, ny-1))
    area.x = array(0, dim = c(nx-1, ny, length(grid.data$Cs_r)))
    area.y = array(0, dim = c(nx, ny-1, length(grid.data$Cs_r)))

    ## Calculate depths
    if (is.na(grid.data$Vtransform) | grid.data$Vtransform == 1) { ## Equation 1
        for (i in 1:length(grid.data$Cs_r)) {
            z0 = (grid.data$s_rho[i] - grid.data$Cs_r[i]) * grid.data$hc + grid.data$Cs_r[i] * (h + zice)
            z[,,i] = z0 + zeta * (1 + z0 / (h + zice)) + zice

            z0 = (grid.data$s_rho[i] - grid.data$Cs_w[i]) * grid.data$hc + grid.data$Cs_w[i] * (h.x + zice.x)
            z.x[,,i] = z0 + zeta.x * (1 + z0 / (h.x + zice.x)) + zice.x

            z0 = (grid.data$s_rho[i] - grid.data$Cs_w[i]) * grid.data$hc + grid.data$Cs_w[i] * (h.y + zice.y)
            z.y[,,i] = z0 + zeta.y * (1 + z0 / (h.y + zice.y)) + zice.y
        }
    } else {
        message(' No valid Vtransform given. Returning zeros!')
    }

    ## Calculate Areas
    if (verbose) { message(' Calculating vertical area panels in x and y...')}
    for (i in length(grid.data$Cs_r):1) {
        if (i == length(grid.data$Cs_r)) {
            for (j in 1:(dim(area.x)[2]-1)) {
                area.x[,j,i] = diff(grid.data$x_v[,j]) * z.x[,j,i]
            }
            j = dim(area.x)[2]
            area.x[,j,i] = diff(grid.data$x_v[,j-1]) * z.x[,j,i]

            for (j in 1:(dim(area.y)[1]-1)) {
                area.y[j,,i] = diff(grid.data$y_u[j,]) * z.y[j,,i]
            }
            j = dim(area.y)[1]
            area.y[j,,i] = diff(grid.data$y_u[j-1,]) * z.y[j,,i]
        } else {
            for (j in 1:(dim(area.x)[2])-1) {
                area.x[,j,i] = diff(grid.data$x_v[,j]) * (z.x[,j,i] - z.x[,j,i+1])
            }
            j = dim(area.x)[2]
            area.x[,j,i] = diff(grid.data$x_v[,j-1]) * (z.x[,j,i] - z.x[,j,i+1])

            for (j in 1:(dim(area.y)[1])-1) {
                area.y[j,,i] = diff(grid.data$y_u[j,]) * (z.y[j,,i] - z.y[j,,i+1])
            }
            j = dim(area.y)[1]
            area.y[j,,i] = diff(grid.data$y_u[j-1,]) * (z.y[j,,i] - z.y[j,,i+1])
        }
    }
    area.x = abs(area.x)
    area.y = abs(area.y)

    if (verbose) { message(' Calculating lateral area panels in x and y...')}
    ## Calculate lateral area of each cell (this is very inefficient but simple, should improve it later!)
    for (i in 1:dim(grid.data$x_u)[1]) {
        for (j in 1:dim(grid.data$y_v)[2]) {
            area.f[i,j] = diff(grid.data$x_v[,j])[i] * diff(grid.data$y_u[i,])[j]
        }
    }

    ## Return
    list(lateral = area.f, x = area.x, y = area.y)
}


#### GRID UTILITIES

#' @export
calc.roms.location = function(lon, lat, grid.data) {

    if (length(lon) != length(lat)) {
        stop('Error, lon and lat length are different!')
    }

    ## Transform to common ranges
    lon = lon %% 360
    grid.data$lon_psi = grid.data$lon_psi %% 360

    locations = data.frame(lon = as.numeric(grid.data$lon_psi), lat = as.numeric(grid.data$lat_psi))
    i = rep(NA, length(lon))
    j = rep(NA, length(lat))

    for (k in 1:length(lon)) {
        a = which.min((locations$lon - lon[k])^2 + (locations$lat - lat[k])^2)
        i[k] = a %% dim(grid.data$lon_psi)[1]
        if (i[k] == 0) { i[k] = dim(grid.data$lon_psi)[1] }
        j[k] = floor((a - i[k]) / dim(grid.data$lon_psi)[1]) + 1
    }

    ## Return
    data.frame(i = i, j = j)
}


#' @export
destagger.grid = function(x) {
    x = as.array(x)
    nx = dim(x)[1]
    ny = dim(x)[2]

    if (length(dim(x)) == 1) {
        x = (x[1:(nx-1)] + x[2:nx]) / 2
    } else if (length(dim(x)) == 2) {
        x = (x[1:(nx-1), 1:(ny-1)] + x[2:nx, 1:(ny-1)] + x[1:(nx-1), 2:ny] + x[2:nx, 2:ny]) / 4
    } else if (length(dim(x)) == 3) {
        x = (x[1:(nx-1), 1:(ny-1),] + x[2:nx, 1:(ny-1),] + x[1:(nx-1), 2:ny,] + x[2:nx, 2:ny,]) / 4
    } else if (length(dim(x)) == 4) {
        x = (x[1:(nx-1), 1:(ny-1),,] + x[2:nx, 1:(ny-1),,] + x[1:(nx-1), 2:ny,,] + x[2:nx, 2:ny,,]) / 4
    } else if (length(dim(x)) == 5) {
        x = (x[1:(nx-1), 1:(ny-1),,,] + x[2:nx, 1:(ny-1),,,] + x[1:(nx-1), 2:ny,,,] + x[2:nx, 2:ny,,,]) / 4
    }

    ## Return x
    x
}

#' @export
destagger.grid.x = function(x) {
    x = as.array(x)
    nx = dim(x)[1]

    if (length(dim(x)) == 1) {
        x = (x[1:(nx-1)] + x[2:nx]) / 2
    } else if (length(dim(x)) == 2) {
        x = (x[1:(nx-1),] + x[2:nx,]) / 2
    } else if (length(dim(x)) == 3) {
        x = (x[1:(nx-1),,] + x[2:nx,,]) / 2
    } else if (length(dim(x)) == 4) {
        x = (x[1:(nx-1),,,] + x[2:nx,,,]) / 2
    } else if (length(dim(x)) == 5) {
        x = (x[1:(nx-1),,,,] + x[2:nx,,,,]) / 2
    }

    ## Return x
    x
}


#' @export
destagger.grid.y = function(x) {
    x = as.array(x)
    nx = dim(x)[2]

    if (length(dim(x)) == 1) {
        nx = dim(x)[1]
        x = (x[1:(nx-1)] + x[2:nx]) / 2
    } else if (length(dim(x)) == 2) {
        x = (x[,1:(nx-1)] + x[,2:nx]) / 2
    } else if (length(dim(x)) == 3) {
        x = (x[,1:(nx-1),] + x[,2:nx,]) / 2
    } else if (length(dim(x)) == 4) {
        x = (x[,1:(nx-1),,] + x[,2:nx,,]) / 2
    } else if (length(dim(x)) == 5) {
        x = (x[,1:(nx-1),,,] + x[,2:nx,,,]) / 2
    }

    ## Return x
    x
}


#' @export
destagger.grid.z = function(x) {
    x = as.array(x)
    nx = dim(x)[3]

    if (length(dim(x)) == 1) {
        nx = dim(x)[1]
        x = (x[1:(nx-1)] + x[2:nx]) / 2
    } else if (length(dim(x)) == 2) {
        nx = dim(x)[2]
        x = (x[,1:(nx-1)] + x[,2:nx]) / 2
    } else if (length(dim(x)) == 3) {
        x = (x[,,1:(nx-1)] + x[,,2:nx]) / 2
    } else if (length(dim(x)) == 4) {
        x = (x[,,1:(nx-1),] + x[,,2:nx,]) / 2
    } else if (length(dim(x)) == 5) {
        x = (x[,,1:(nx-1),,] + x[,,2:nx,,]) / 2
    }

    ## Return x
    x
}





#####################
#### OLD STUFF
#####################

#' @title Get ROMS Time old
#' @export
#' @author Thomas Bryce Kelly
#' @import ncdf4
get.roms.time.old = function(path, convert = TRUE) {

    ## Open the nc file
    roms = ncdf4::nc_open(path, write=FALSE)

    nhist = ncdf4::ncvar_get(roms, varid = 'nHIS') # dt * nHist = time between records
    dstart = ncdf4::ncvar_get(roms, varid = 'dstart') # days since 1900-01-01 00:00:00
    dt = ncdf4::ncvar_get(roms, varid = 'dt') # Sec
    ntimes = ncdf4::ncvar_get(roms, varid = 'ntimes') # number of time frames

    ## Close the nc file
    ncdf4::nc_close(roms)

    ## Calculate the times (seconds from 1900-01-01)
    roms.times = seq(from = dstart, by = dt * nhist / 86400, length.out = ntimes/nhist) * 86400

    if (convert) {
        roms.times = conv.time.roms(roms.times, tz='GMT')
    }
    roms.times
}



#' @title Get ROMS z-levels
#' @author Thomas Bryce Kelly
calc.roms.depths.old = function(h, hc, theta, b = 0.6, N) {
    z = array(NA, dim = c(dim(h), N))
    ds = 1 / N
    s = seq(-1 + ds/2, -ds/2, by = ds)
    A = sinh(theta * s) / sinh(theta)
    B = (tanh(theta * (s + 0.5)) - tanh(theta * 0.5)) / (2 * tanh(theta * 0.5))
    C = (1 - b) * A + b * B

    for (i in 1:dim(h)[1]) {
        for (j in 1:dim(h)[2]) {
            z[i,j,] = hc * s + (h[i,j] - hc) * C
        }
    }
    z
}

#' @export
load.roms.old = function(path, verbose = T) {
    if (verbose) { message(Sys.time(), ': Attempting to load ROMS file ', path, appendLF = F)}
    file = ncdf4::nc_open(path, write=FALSE)
    if (verbose) { message(' success.')}

    vars = names(file$var)

    #### Load variables
    ## Vertical Grid
    if ('h' %in% vars) {
        if (verbose) { message(' Loading h...')}
        h = ncdf4::ncvar_get(file, 'h')
    } else {
        stop('Error, no h variable found.')
    }

    if ('hc' %in% vars) {
        if (verbose) { message(' Loading hc...')}
        hc = ncdf4::ncvar_get(file, 'hc')
    } else {
        stop('Error, no hc variable found.')
    }

    if ('theta_s' %in% vars) {
        if (verbose) { message(' Loading theta.s...')}
        theta.s = ncdf4::ncvar_get(file, 'theta_s')
    } else {
        stop('Error, no theta_s variable found.')
    }

    if ('theta_b' %in% vars) {
        if (verbose) { message(' Loading theta.b...')}
        theta.b = ncdf4::ncvar_get(file, 'theta_b')
    } else {
        message('Warning, no theta_b variable found, setting to zero.')
        theta.b = 0
    }

    h = (h[2:dim(h)[1],] + h[1:(dim(h)[1] - 1),]) / 2 ## interpolate to center of each grid cell
    h = (h[,2:dim(h)[2]] + h[,1:(dim(h)[2] - 1)]) / 2

    ## Horizontal
    if ('lat_psi' %in% vars) {
        if (verbose) { message(' Loading lat_psi...')}
        lat = ncdf4::ncvar_get(file, 'lat_psi')

        if (verbose) { message(' Loading lon_psi...')}
        lon = ncdf4::ncvar_get(file, 'lon_psi')
    } else if ('x_psi' %in% vars) {
        if (verbose) { message(' Loading x_psi...')}
        lon = ncdf4::ncvar_get(file, 'x_psi')
        if (verbose) { message(' Loading y_psi...')}
        lat = ncdf4::ncvar_get(file, 'y_psi')
    } else {
        message('No positional coordiantes found!')
    }

    grid = list(lon = lon, lat = lat)
    if (verbose) { message(' Formed grid ', length(grid$lon), 'x', length(grid$lat))}

    ## Advect
    if ('u' %in% vars) {
        if (verbose) { message(' Loading u... ', appendLF = F); a = Sys.time()}
        u = ncdf4::ncvar_get(file, 'u')
        if (verbose) { message('(',Sys.time() - a, ')')}
        u = (u[,2:dim(u)[2],,] + u[,1:(dim(u)[2] - 1),,]) / 2
    } else if ('Huon' %in% vars) {
        if (verbose) { message(' Loading Huon... ', appendLF = F); a = Sys.time()}
        u = ncdf4::ncvar_get(file, 'Huon')
        if (verbose) { message('(',Sys.time() - a, ')')}
        u = (u[,2:dim(u)[2],,] + u[,1:(dim(u)[2] - 1),,]) / 2
    } else {
        stop('Error, no u velocity found.')
    }

    if ('v' %in% vars) {
        if (verbose) { message(' Loading v... ', appendLF = F); a = Sys.time()}
        v = ncdf4::ncvar_get(file, 'v')
        if (verbose) { message('(',Sys.time() - a, ')')}
        v = (v[2:dim(v)[1],,,] + v[1:(dim(v)[1] - 1),,,]) / 2
    } else if ('Hvom' %in% vars) {
        if (verbose) { message(' Loading Hvom... ', appendLF = F); a = Sys.time()}
        v = ncdf4::ncvar_get(file, 'Hvom')
        if (verbose) { message('(',Sys.time() - a, ')')}
        v = (v[2:dim(v)[1],,,] + v[1:(dim(v)[1] - 1),,,]) / 2
    } else {
        stop('Error, no v velocity found.')
    }

    if ('w' %in% vars) {
        if (verbose) { message(' Loading w... ', appendLF = F); a = Sys.time()}
        w = ncdf4::ncvar_get(file, 'w')
        if (verbose) { message('(',Sys.time() - a, ')')}
        w = (w[,,2:dim(w)[3],] + w[,,1:(dim(w)[3] - 1),]) / 2 # depth averaged w
        w = (w[2:dim(w)[1],,,] + w[1:(dim(w)[1] - 1),,,]) / 2
        w = (w[,2:dim(w)[2],,] + w[,1:(dim(w)[2] - 1),,]) / 2
    } else {
        stop('Error, no w velocity found.')
    }

    if (verbose) { message(' Unifying grid layout...')}
    dims = c(dim(u)[1] * dim(u)[2], dim(u)[3], dim(u)[4])
    dim(u) = dims
    dim(v) = dims
    dim(w) = dims
    dim(lon) = dims[1]
    dim(lat) = dims[1]
    dim(h) = dims[1]
    N = dims[2]

    ## Diffusivity
    if ('AKt' %in% vars) {
        if (verbose) { message(' Loading AKt...')}
        AKt = ncdf4::ncvar_get(file, 'AKt')
    } else {
        AKt = NULL
    }

    ## Fields
    if ('temp' %in% vars) {
        if (verbose) { message(' Loading pot temp...')}
        temp = ncdf4::ncvar_get(file, 'temp')
    } else {
        temp = NULL
    }

    if ('rho' %in% vars) {
        if (verbose) { message(' Loading rho...')}
        rho = ncvar_get(file, 'rho')
    } else {
        rho = NULL
    }
    if ('salt' %in% vars) {
        if (verbose) { message(' Loading salinity...')}
        salt = ncvar_get(file, 'salt')
    } else {
        salt = NULL
    }

    if (verbose) { message(' Adjusting tracer grid...')}
    if (!is.null(temp)) {
        temp = (temp[2:dim(temp)[1],,,] + temp[1:(dim(temp)[1] - 1),,,]) / 2
        temp = (temp[,2:dim(temp)[2],,] + temp[,1:(dim(temp)[2] - 1),,]) / 2
        dim(temp) = dims
    }

    if (!is.null(rho)) {
        rho = (rho[2:dim(rho)[1],,,] + rho[1:(dim(rho)[1] - 1),,,]) / 2
        rho = (rho[,2:dim(rho)[2],,] + rho[,1:(dim(rho)[2] - 1),,]) / 2
        dim(rho) = dims
    }

    if (!is.null(salt)) {
        salt = (salt[2:dim(salt)[1],,,] + salt[1:(dim(salt)[1] - 1),,,]) / 2
        salt = (salt[,2:dim(salt)[2],,] + salt[,1:(dim(salt)[2] - 1),,]) / 2
        dim(salt) = dims
    }

    if (!is.null(AKt)) {
        AKt = (AKt[2:dim(AKt)[1],,,] + AKt[1:(dim(AKt)[1] - 1),,,]) / 2
        AKt = (AKt[,2:dim(AKt)[2],,] + AKt[,1:(dim(AKt)[2] - 1),,]) / 2
        AKt = (AKt[,,2:dim(AKt)[3],] + AKt[,,1:(dim(AKt)[3] - 1),]) / 2
        dim(AKt) = dims
    }

    #### Calculated
    z = get.zlevel.roms(h, hc, theta.s, N=N)
    time = get.roms.time(path, TRUE)
    meta = list(path = path, theta.b = theta.b, theta.s = theta.s, hc = hc, N = N, dims = dims)

    #### Return object
    list(lat = lat,
         lon = lon,
         grid = grid,
         depth = function(x) {-get.zlevel.roms(x, hc, theta.s, N = N)},
         z = z,
         h = h,
         u = u,
         v = v,
         w = w,
         AKt = AKt,
         T = temp,
         rho = rho,
         S = salt,
         time = time,
         meta = meta)
}



#' @title Load LTRANS Data
#' @author Thomas Bryce Kelly
#' @export
#' @import ncdf4
load.ltrans = function(path, roms.path, backward = TRUE) {
    cycle = ncdf4::nc_open(path, write=FALSE)
    lat = ncdf4::ncvar_get(cycle, 'lat')
    lon = ncdf4::ncvar_get(cycle, 'lon')
    S = ncdf4::ncvar_get(cycle, 'salinity')
    T = ncdf4::ncvar_get(cycle, 'temperature')
    age = ncdf4::ncvar_get(cycle, 'age')
    dob = ncdf4::ncvar_get(cycle, 'dob')
    depth = ncdf4::ncvar_get(cycle, 'depth')
    model.time = ncdf4::ncvar_get(cycle, 'model_time')
    ncdf4::nc_close(cycle)

    rt = get.roms.time(roms.path, FALSE)

    if (backward) {
        model.time = max(rt) - model.time
        dob = max(rt) - dob
    }

    ## Convert to POSIX
    model.time = conv.time.roms(model.time, tz='GMT')
    dob = conv.time.roms(dob, tz='GMT')
    rt = conv.time.roms(rt, tz='GMT')

    ## Return
    list(file = path, lat = lat, lon = lon, sal = S, temp = T, n.times = length(model.time),
                 age = age, dob = dob, depth = depth, time = model.time, n.particles = length(lon[,1]),
                 roms.times = rt)
}