I started to tackle this in #641 because I wanted to have the derivative functions (e.g. drdn_dsurface) call the other surface methods (e.g. calc_rfl) directly. This is tricky however, because we have to track the potential (mis)match between the RT and Surface wavelength grids.

Currently, the up-sampling is done within forward.py, which makes a lot of sense, but there is a lot of redundancy, and in the case of the calc_rfl function in surface_glint_model.py, there doesn't look like there is a check to avoid len(x_surface_rfl) != len(L_down_dir). Do these situations exist? If so, the upsampling in forward.py is performed after calc_rfl is called, and the addition of the sky glint term would raise an error.

        rho_dir_dir, rho_dif_dir = self.surface.calc_rfl(
            x_surface, geom, L_down_dir, L_down_dif
        )
        rho_dir_dir_hi = self.upsample(self.surface.wl, rho_dir_dir)
        rho_dif_dir_hi = self.upsample(self.surface.wl, rho_dif_dir)

This is just a thought to revisit the upsampling throughout the forward.py, surface, RT, and instrument.