pam_fmado_ll <- function (x, k, ll) {

    # x - Design matrix, typically colums are stations, rows are time 
    # k - Number of clusters
    # ll - two column matrix of lat long points (or preferably projected) with N rows
    
    N = ncol(x) # number of stations 
    T = nrow(x) # number of time points

    # compute the F-madogram distance
    V = array(NaN, dim = c(T, N))
    for (p in 1:N) {
        x.vec = as.vector(x[, p])
        Femp = ecdf(x.vec)(x.vec)
        V[, p] = Femp
    }
    DD = dist(t(V), method = "manhattan", diag = TRUE, upper = TRUE)/(2 * T)

    # weight by physical distance
    DDll = dist(ll,method='manhattan')
    DDw = as.matrix(DD) + t(t(as.matrix(DDll))/apply(as.matrix(DDll),2,max))*max(as.matrix(DD))

    # do the clustering
    output = pam(DDw, k, diss = TRUE, medoids = NULL)
    return(output)
}