docs: def interface_temperature(Ts, Tatm, **kwargs):
      '''Compute temperature at model layer interfaces.'''
      #  Actually it's not clear to me how the RRTM code uses these values
      lev = Tatm.domain.axes['lev'].points
      lev_bounds = Tatm.domain.axes['lev'].bounds
      #  Interpolate to layer interfaces
      f = interp1d(lev, Tatm, axis=-1)  # interpolation function
      Tinterp = f(lev_bounds[1:-1])
      #  add TOA value, Assume surface temperature at bottom boundary
      Ttoa = Tatm[...,0]
      Tinterp = np.concatenate((Ttoa[..., np.newaxis], Tinterp, Ts), axis=-1)
      return Tinterp