Returns layer depth at lat / lon (degrees)
    where lat/lon may be arrays (of equal size).
    Depths are returned in metres.

def layer_depth( lat, lon, layerID="LID-BOTTOM"):
    """Returns layer depth at lat / lon (degrees)
    where lat/lon may be arrays (of equal size).
    Depths are returned in metres.
    """

    ## Must wrap longitude from 0 to 360 ...

    lon1 = np.array(lon)%360.0
    lat1 = np.array(lat)

    # ## Must wrap longitude from -180 to 180 ...
    #
    # lon1[np.where(lon1 > 180.0)] = 360.0 - lon1[np.where(lon1 > 180.0)]
    #
    data, err = _interpolator.interpolate( np.radians(lon1), np.radians(lat1),
                                      _litho_data[l1_layer_decode[layerID], l1_data_decode["DEPTH"]], order=1 )

    return data