def divergence(U, V, pm, pn):
    """Divergence from C-grid velocity

    U : U-component, velocity [m/s], shape[-2:] = (jmax, imax+1)
    V : V-component, velocity [m/s], shape[-2:] = (jmax+1, imax)
    pm : invers metric term X-direction [1/m], shape = (jmax+2, imax+2)
    pn : invers metric term Y-direction [1/m], shape = (jmax+2, imax+2)

    Result:
    div : divergence [1/s],  shape[-2:] = (jmax, imax), [1/s]

    The divergence for an orthogonal coordinate system is given by::

      div = pm * pn * (d/dx(U/pn) + d/dy(V/pm))

    """

    A = 2.0 * U / (pn[1:-1, :-1] + pn[1:-1, 1:])  # U/pn
    B = 2.0 * V / (pm[:-1, 1:-1] + pm[1:, 1:-1])  # V/pm
    pmn = pm[1:-1, 1:-1] * pn[1:-1, 1:-1]

    div = pmn * (A[..., :, 1:] - A[..., :, :-1] + B[..., 1:, :] - B[..., :-1, :])

    return div