Linear interpolation issue - Meta.top

[231Harriet]

Hi,

I am interested in the interpolation method used in the meta.top function. In the Read et al., (2011) paper it states that the method is designed to; "interpolate between i and i+1 to yield the approximate depth of the base of the mixed layer". That is, the method interpolates between the depth where the threshold is first met and the consecutive depth. However, I find in the code the interpolation is actually between the depth before the threshold was met and the depth of the depth of the thermocline. Is this possibly a mistake?

This is the specific line of code:

metaTop_depth = stats::approx(drho_dz[i:thermo_index], sortDepth[i:thermo_index], slope)

which could be changed to:

metaTop_depth = stats::approx(drho_dz[i:i+1], sortDepth[i:i+1], slope)

This is the original full function:

function (wtr, depths, slope = 0.1, seasonal = TRUE, mixed.cutoff = 1) 
{
    if (any(is.na(wtr))) {
        return(rep(NaN, 2))
    }
    if (length(wtr) < 3) {
        return(c(max(depths), max(depths)))
    }
    depths = sort.int(depths, index.return = TRUE)
    wtr = wtr[depths$ix]
    depths = depths$x
    thermoD = thermo.depth(wtr, depths, seasonal = seasonal, 
        mixed.cutoff = mixed.cutoff)
    if (is.na(thermoD)) {
        return(c(NaN, NaN))
    }
    rhoVar = water.density(wtr)
    dRhoPerc = 0.15
    numDepths = length(depths)
    drho_dz = vector(mode = "double", length = numDepths - 1)
    for (i in 1:(numDepths - 1)) {
        drho_dz[i] = (rhoVar[i + 1] - rhoVar[i])/(depths[i + 
            1] - depths[i])
    }
    metaBot_depth = depths[numDepths]
    metaTop_depth = depths[1]
    Tdepth = rep(NaN, numDepths - 1)
    for (i in 1:(numDepths - 1)) {
        Tdepth[i] = mean(depths[i:(i + 1)])
    }
    tmp = sort.int(unique(c(Tdepth, thermoD)), index.return = TRUE)
    sortDepth = tmp$x
    sortInd = tmp$ix
    numDepths = length(sortDepth)
    drho_dz = stats::approx(Tdepth, drho_dz, sortDepth)
    drho_dz = drho_dz$y
    thermo_index = which(sortDepth == thermoD)
    for (i in thermo_index:numDepths) {
        if (!is.na(drho_dz[i]) && drho_dz[i] < slope) {
            metaBot_depth = sortDepth[i]
            break
        }
    }
    if (i - thermo_index >= 1 && (!is.na(drho_dz[thermo_index]) && 
        drho_dz[thermo_index] > slope)) {
        metaBot_depth = stats::approx(drho_dz[thermo_index:i], 
            sortDepth[thermo_index:i], slope)
        metaBot_depth = metaBot_depth$y
    }
    if (is.na(metaBot_depth)) {
        metaBot_depth = depths[numDepths]
    }
    for (i in seq(thermo_index, 1)) {
        if (!is.na(drho_dz[i]) && drho_dz[i] < slope) {
            metaTop_depth = sortDepth[i]
            break
        }
    }
    if (thermo_index - i >= 1 && (!is.na(drho_dz[thermo_index]) && 
        drho_dz[thermo_index] > slope)) {
        metaTop_depth = stats::approx(drho_dz[i:thermo_index], 
            sortDepth[i:thermo_index], slope)
        metaTop_depth = metaTop_depth$y
    }
    if (is.na(metaTop_depth)) {
        metaTop_depth = depths[i]
    }
    return(c(metaTop_depth, metaBot_depth))
}

Thanks,
Harriet