If $n$ is even, $f$ has to be the sum $\sum_{i=1}^{n}{|p_i-c|}$ or a monotonically increasing function of this sum.
Your requirement: $f(x_1,...,x_n)\le f(y_1,...,y_n)$ whenever $x_i=|p_i-c|$, $y_i=|p_i-d|$, and $c$ is median of $p$, written as $c=m(p)$. The vector $x$ represents distances from median location and vector $y$ represents distances from another location.
Note that $x_i$ and $y_i$ cannot be arbitrary vectors. Suppose $\delta=|c-d|$. If $\delta<\mathrm{min}_{i\neq j}{\{x_i-x_j\}}$, then $x_i=y_i+\delta$ for at least $\lceil n/2 \rceil$ coordinates and $x_i=y_i-\delta$ for remaining coordinates. If $\delta>\mathrm{min}_{i\neq j}{\{x_i-x_j\}}$, then $x_i=y_i+\delta$ for at least $\lceil n/2 \rceil$ coordinates and $x_i \ge y_i-\delta$ for remaining coordinates.
So $f$ has to be a function which increases when all its arguments are changed by a small number $\delta>0$ and the number of arguments increased exceeds the number of arguments decreased. In the following, I assume $f$ is differentiable in all arguments.
Consider $n$ even. At $x$, all first partial derivatives of $f$ must be equal: $\frac{\partial f}{\partial x_i} = \frac{\partial f}{\partial x_j} \forall i,j$ otherwise we could increase $n/2$ arguments with lowest partial derivatives and decrease the remaining arguments and this will lead to a lower value of $f.$ Moreover, $x$ can be any set of coordinates (median doesn't have to be one of the points). So all partial derivatives of $f$ are equal everywhere in domain of $f$. It is easy to show then that increasing $x_i$ and decreasing $x_j$ by the same amount leaves $f$ unchanged. So $f$ depends on $x_i+x_j$ rather than $x_i$ and $x_j$ separately. Since, $i$ and $j$ are arbitrary, $f$ is a weakly increasing function of $\sum_{i=1}^{n}x_i$. Note that $f$ doesn't have to be an affine function of the sum.
If $n$ is odd, $k \ge 1$ coordinates must be zero in $x$. Then, a small change from $x$ to $y$ increases these $k$ coordinates by $\delta$. Of the remaining coordinates, some (at least $l=\frac{n+1}{2}-k$ and at most $\frac{n-1}{2}$) will increase by $\delta$ while the remaining will decrease by $\delta$. Since more than half of the coordinates increase, function $f$ will increase if all partial derivatives are equal and positive. That is, $f$ can again be an increasing function of $\sum_{i=1}^{n}x_i$. However, this is not the only possibility. There may be other solutions which satisfy the following necessary condition:
If $k$ coordinates of $x$ are $0$, reorder the coordinates and accordingly the arguments of $f$ so that $x_i=0$ for $1 \le i \le k$, $x_i>0$ for $k < i \le n$ and $\frac{\partial f}{\partial x_i} \le \frac{\partial f}{\partial x_j} \forall k<i<j$. Then, $\sum_{i=1}^{k+l}{\frac{\partial f}{\partial x_i}} \ge \sum_{i=k+l+1}^{n}{\frac{\partial f}{\partial x_i}}$.
The above necessary condition ensures that $f$ has a local minimum at $x$ but need not guarantee a global minimum. The functions $\prod_{i=1}^{n}{(x_i)}$ and $\min_{1\le i \le n}{(x_i)}$ satisfy the necessary condition.
Finally, note that the minimum and product functions do not work for the case of even $n$. For example, try $p=(1,3,11,13,15,39)$.
The problem for odd $n$ can be characterized as follows. Define an ordered pair of vectors $(x,y)$ as comparable if $x$ represents distances from median in a network and $y$ represents distances from another point in the same network. Let $\Omega$ represent the set of comparable points. Then, $f$ satisfies
Condition A: $f(x) \le f(y) \quad \forall (x,y) \in \Omega$.
While there are no direct restrictions on relative values of $f$ for vectors $x$ and $y$ if $(x,y) \not\in \Omega$, repeated application of Condition A may yield such restrictions. For example, if $(x,z) \in \Omega$ and $(z,y) \in \Omega$ then $f(x) \le f(y)$. Characterizing desired functions requires understanding the structure of $\Omega$.
$(x,y) \in \Omega$ if and only if at least one coordinate of $x$ is zero and there exists $\delta \in \mathbb{R}$ such that there are $\lceil n \rceil$ values of coordinate $i$ such that $y_i=|x_i+\delta|$ and $\lceil n \rceil$ values of coordinate $i$ such that $y_i=|x_i-\delta|$. The trivial solutions you mentioned equate $f$ to its lowest value at all $x$ with a coordinate of $0$. But this means that if $(x,y) \in \Omega$ for some $y$, then $f(x) \le f(z) \forall z$ even when $(x,z) \not\in \Omega$. Increasing functions of sum of coordinates exploit the property that if $(x,y) \in \Omega$ then $\sum x_i < \sum y_i$. However, note that $\sum x_i < \sum y_i \not\Rightarrow (x,y) \in \Omega$.
One approach towards finding other functions meeting Condition A would be to characterize the set of all pairs $(x,y)$ of median-distance vectors (at least one component of each of $x$ and $y$ is zero) at which the realtive values of $f$ are restricted. That is, all pairs $(x,y)$ such that there exist $x=z_0,z_1,z_2,...,z_{m-1},z_m=y$ with $(z_{i-1},z_i) \in \Omega$. This step seems difficult. After this step, choose function values at all median-distance vectors satisfying these restrictions. Finally, fill in values at other vectors $y$ that are not distances from medians to exceed $f(x)$ for all $x$ such that $(x,y) \in \Omega$.
Since only relative function values are restricted only for $(x,y) \in \Omega$, it seems unlikely that $f$ can be characterized by properties which depend solely on how values of $f$ change in a neighborhood.
