code: def get_layer_nr(inp, depth):
      r
      zinp = inp[2]
  
      #  depth = [-infty : last interface]; create additional depth-array
      # pdepth = [fist interface : +infty]
      pdepth = np.concatenate((depth[1:], np.array([np.infty])))
  
      # Broadcast arrays
      b_zinp = np.atleast_1d(zinp)[:, None]
  
      # Get layers
      linp = np.where((depth[None, :] < b_zinp)*(pdepth[None, :] >= b_zinp))[1]
  
      # Return; squeeze in case of only one inp-depth
      return np.squeeze(linp), zinp
docs: r"""Get number of layer in which inp resides.
  
      Note:
      If zinp is on a layer interface, the layer above the interface is chosen.
  
      This check-function is called from one of the modelling routines in
      :mod:`model`.  Consult these modelling routines for a detailed description
      of the input parameters.
  
      Parameters
      ----------
      inp : list of floats or arrays
          Dipole coordinates (m)
  
      depth : array
          Depths of layer interfaces.
  
  
      Returns
      -------
      linp : int or array_like of int
          Layer number(s) in which inp resides (plural only if bipole).
  
      zinp : float or array
          inp[2] (depths).