There are two ways to consider the solvation, through an implicit and explicit method. The first is based on considering the solvation through a conducted-like polarizable method, such as CPCM or universal model (SMD). I think SMD is the best for electronic states. The second is based on adding the solvent molecules around the compound and performing a multilevel QM/MM, where the compound is treated at the quantum-mechanics level and the solvent is treated under molecular mechanics, or even QM/QM with lower level QM for the solvent, such as GFNn-xTB.
In my conception, I think it depends on the information you are trying to extract, but the first is less computationally expansive and considers only the refractive index and the dielectric constant of the solvent. You can use the Gaussian software or Orca to calculate it.
Got it. It all depends on the system you are modeling because if you are dealing with an organic molecule, the B3LYP/6-31+g(d,p) is ok at best. However, for an inorganic system such as a coordination compound, you have to use a larger basis set. Also, I don't recommend the use of the PCM solvation model, CPCM or SMD are better models.
Here is an example of a robust procedure:
Optimization: B3LYP D4/Def2-SVP for organic molecules, where D4 stands for the dispersion correction, and B3LYP D4/Def2-TZVP for coordination compounds.
Single point: PBE0 D4/Def2-SVP or PBE0 D4/Def2-TZVP.
If you are dealing with excited states, I recommend the use of range-separates hybrid functionals, such as CAM-B3LYP or wB97X.