How is Meissner effect explained by BCS theory?
The bottom line is the spontaneous symmetry breakdown from global $U(1)$ to $\mathbb{Z}_2$ and the concomitant rigidity of the omnipresent coherent phase down to which the system breaks. However, both the microscopic action and the BCS ground state (3) of a superconductor possess local $U(1)$ gauge symmetry.
By rigidity, I mean something reminiscent of a restoring force felt when one tries bending or distorting a solid stick, which, fundamentally, originates from the translational symmetry breaking in a typical crystalline solid. Note that this global $U(1)$ symmetry breakdown does not contradict Elitzur's theorem which forbids spontaneous breakdown of local gauge symmetry. The value of this phase in a superconductivity ground state is not observable (maximal uncertainty) since one has to integrate over $\phi\in[0,2\pi)$ so as to recover particle conservation. Nonetheless, its variation in spacetime, hence rigidity, is not only physically observable but also crucial, especially in Meissner effect and Josephson effect.
The form of BCS ground state (3) (given at the end) says, all Cooper pairs of different momenta $\vec{k}$s share exactly the same in-pair relative phase $\phi\equiv \mathrm{arg}(v_{\vec{k}}/u_{\vec{k}})$, namely the aforesaid omnipresent coherent phase is almost of the same value inside the whole superconductor. So does the gap function $\Delta$ in the mean-field BCS interaction term $\Delta c^{\dagger}\left(x\right)c^{\dagger}\left(x\right)+\textrm{h.c.}$. This also manifests the symmetry breaking down to $\mathbb{Z}_2$ for that only $\{0,\pi\}$ phase transformation of $\{c,c^\dagger\}$ renders the term invariant. (See also this good answer, owing to a discussion with the author of which, I managed to correct the inaccuracies in my answer.)
In the superconducting phase, one can choose the unitary gauge to make the Goldstone field $\phi(x)=0$ everywhere and certainly arrive at some new $A'_\mu$ (physically unnecessary, however, having the merit of vanishing Goldstone modes and no Faddeev-Popov ghost). Or one can instead replace the massless gauge field $A_\mu$ by a newly defined non-gauge vector field $-\frac{e^*}{c}A'_\mu\equiv \hbar\partial_\mu\phi-\frac{e^*}{c}A_\mu$ while preserving the total three degrees of freedom. Anyway, the new $A'_\mu$ is manifestly massive in the effective Lagrangian. And it looks as if photons coupled to a superconductor undergo kinda explicit gauge symmetry breaking and acquire mass via this Abelian Higgs mechanism, which restrains the long-range electromagnetic field to sort of exponentially decaying Yukawa potential. This is nothing but the Meissner effect clarified below Eq.(1).
And consequently, a macroscopic coherent quantum state with phase rigidity, just as (3) does, is constructed. Equivalently speaking, this is a phenomenon where the quantum mechanical phase reaches macroscopic dimensions, which is somewhat natural for Bosons (e.g. Bose-Einstein condensation and Bosonic superfluidity of ${}^4\mathrm{He}$), and is astonishingly also achieved via the formation of Fermions' Cooper pairs.
If you need more calculations
Rigidity dictates Meissner effect. As is well known, in the presence of electromagnetic field (no matter whether it is a penetrating one or a conventional one), the wavefunction gains an Aharonov-Bohm $U(1)$ phase $\mathrm{e}^{\mathrm{i}e\chi(x)/\hbar}$, wherein nonintegrable phase $\chi=\int{A_\mu\mathrm{d}x^\mu}$ might, in general, be path-dependent. What if this system possesses some rigidity of this distribution of twisting angle $\chi(x)$? In analogy to the twisting or distortion in a solid body, a macroscopic term $\int{\frac{1}{2}\kappa(\nabla\chi)^2\mathrm{d}V}$ arises in the free energy of this system. More detailed analysis in BCS theory indeed gives you a free energy increase (c.f. last section of this answer) $$ \Delta G=e^2\frac{\rho}{2m}\int{A^2}\mathrm{d}V, $$ wherein electron density $\rho=\langle\psi^*(x)\psi(x)\rangle$ (spin degrees of freedom neglected for the nonce). And as a result of this, we get one of the London's equations $$ \vec{j}_d=-\frac{\delta\Delta G}{\delta\vec{A}}=-e^2\frac{\rho}{m}\vec{A},\tag{1} $$ which is the famous $j\propto A$ relation. Combined with Maxwell Eq. $\vec{j}=\nabla\times\vec{H}$ (provided $\vec{A}$ has no temporal dependence), you can easily obtain the exponential decay of $\vec{A}$ or $\vec{B}$ inside the superconductor, that is to say, Meissner effect is mandatorily required by this rigidity. In a nutshell, superconductivity serves as a mechanism that resists the generation of Aharonov-Bohm phase due to penetrated electromagnetic field.
Why does the diamagnetic current persist? In quantum mechanics, electrical current equals to the paramagnetic current subtracted by the diamagnetic current $$ \vec{j}=\vec{j}_p+\vec{j}_d\equiv [\frac{1}{2m}(\psi^*\hat{\vec{p}}\psi-\psi\hat{\vec{p}}\psi^*)]+[-\frac{q}{m}\vec{A}\psi^*\psi].\tag{2} $$
- In a normal state, presence of $\vec{A}$ also increases the free energy, however, in a relative banal way, that is, $\Delta G=\frac{1}{2}\chi\int{(\nabla\times\vec{A})^2\mathrm{d}V}$. Together with Maxwell equations, only a small Landau diamagnetic current is retained $\vec{j}=-\frac{\delta\Delta G}{\delta\vec{A}}=\nabla\times\vec{M}$, where $\vec{M}$ is the local magnetization. This is because the $\vec{j}_p$ and $\vec{j}_d$ in Eq. (2) cancel each other, as is straightforward to check once you notice that $\psi$ contains the aforementioned Aharonov-Bohm phase $\mathrm{e}^{\mathrm{i}e\chi(x)/\hbar}$.
- On the other hand, there's no such cancellation in a superconducting phase. Paramagnetic current $\vec{j}_p$ obviously contains some spatial derivative of the phase in $\psi(x)$, i.e., kinda strain of the wavefunction. However, such twisting is certainly not energetically favoured since the previously discussed rigidity. With hindsight, you might even think it in a sloppy way -- rigidity repels $\vec{A}$ out, no $\vec{A}$ no Aharonov-Bohm phase in $\psi$, $\vec{j}_p=0$ consequently. Anyway, diamagnetic current $\vec{j}_d$ in (1) perfectly persists in the end. (Partly cancelled by nonzero $\vec{j}_p$ when $0<T<T_c$ actually.)
BCS theory provides the microscopic mechanism that yields this rigidity.
- From a field theoretic point of view of BCS theory, we can introduce auxiliary Bosonic fields $\Delta\equiv |\Delta(x)|\mathrm{e}^{\mathrm{i}\theta(x)},\bar{\Delta},\varphi$ so as to perform Stratonovich-Hubbard transformation onto the BCS action $S$. Afterwards, part of the action reads $\frac{1}{2m}\sum_\sigma{\int{(\nabla\theta(x))^2\bar{\psi}_\sigma(x)\psi_\sigma(x)\mathrm{d}V}}$, wherein $\psi$ is the original fermionic field, which conspicuously manifests rigidity of phase $\theta$. Indeed, $\Delta$ corresponds to the superconducting gap or the order parameter. Further, after tedious manipulations and approximations to construct an effective low-energy theory of $\varphi$ and $\theta$, we can directly calculate to show that the paramagnetic current density correlation function $\langle j_p^\alpha(x)j_p^\beta(0)\rangle$ vanishes when $T=0$, which undoubtedly endorses our previous discussion in section 2. You can even see that this owes to the existence of the gap $\Delta$.
- To connect the aforesaid phase $\theta$ of $\Delta$ with the phase of wavefunction, we might turn to the BCS ground state $$ \vert\Psi_{\mathrm{BCS}}(\phi)\rangle=\prod_{\vec{k}}{(|u_{\vec{k}}|+|v_{\vec{k}}|\mathrm{e}^{\mathrm{i}\phi}c_{\vec{k}\uparrow}^\dagger c_{-\vec{k}\downarrow}^\dagger)\vert 0\rangle}.\tag{3} $$ In either BCS's original variational calculation or Bogoliubov transformation approach, this relative phase $\phi\equiv \mathrm{arg}(v_{\vec{k}}/u_{\vec{k}})$ is always directly related to gap $\Delta$ because of $\Delta_\vec{k}^*v_\vec{k}/u_\vec{k}\in \mathbb{R}$. At this stage, we can once again say that the wave function becomes solid and no strain occurs, therefore $\vec{j}_p$ does not emerge.
In all theories (London, Ginzburg-Landau, Bardeen-Cooper-Schrieffer, Bogoliubov-Gennes) the Meissner effect is accounted for as the constitutive relation $j\propto A$ : the current is proportional to the vector-potential. (The proportionality factor depends on the system of units, so let us forget about that here.)
For the London's phenomenological theory, this relation is chosen in order to verify the Meissner relation, and is called the London's (constitutive) relation in his honour. In semi-phenomenological theory à la Ginzburg-Landau, that's the gauge structure in addition to the phase-transition which forces the proportionality between $j$ and $A$.
In microscopic theory, this relation is a bit fake, since the Schrödinger equation already has this structure. Recall the current is $j\propto \Im\left(\Psi \nabla \Psi^{\ast}\right)+A\left|\Psi\right|^{2}$ in (first quantised) quantum mechanics. So what you have to show is that the paramagnetic contribution ($j_{P}\propto\Im\left(\Psi \nabla \Psi^{\ast}\right)$) disappears in a superconductor, and only the diamagnetic contribution $j_{D}\propto A\left|\Psi\right|^{2}$ survives. Note that in the Ginzburg-Landau theory the paramagnetic term disappears by a "gauge-choice". In modern terminology this "gauge-choice" is called an Anderson-Higgs mechanism. From microscopic theory you can prove that the paramagnetic contribution vanishes.
All the calculations are crystal-clear in the original paper
J. Bardeen, L. N. Cooper, and J. R. Schrieffer, Theory of Superconductivity Phys. Rev. 108, 1175 (1957).
which is accessible for free from the editor's website. See especially section V. Electrodynamics properties.
The proportionality factor between $j$ and $A$ has something to do with the London's penetration length $\Lambda$: $j\propto A/\Lambda$ (valid for some superconductors only, please find the relevant details in the BCS's paper). This penetration length can be directly measured, and is proportional to the amplitude of the gap, the superconducting order parameter. That's the correct prediction of the generic trends of the penetration length which served as (one of) the basis for the verification of the BCS theory.
Actually, what BCS proved is that $\lim_{q\rightarrow 0}j\propto A$ since they considered the low-energy sector. From the expression of the first-quantised current, a Fourier transform shows that $\lim_{q\rightarrow 0}j\propto A$ in general, and so it seems that the Meissner effect is a generic property of quantum systems. So the question: do we really need BCS to understand superconductivity? The answer is clearly yes, and the reason follows. In quantum field theory, only bosons have such a Schrödinger-like expression for the current. So in short the Meissner effect is a natural effect for bosons. Conclusion: what is really fundamental for superconductivity is not the $j\propto A$ relation (since this is always valid for any bosons systems in the limit $q\rightarrow 0$), but to understand how a fermionic liquid behaves as a bosonic one !?! So in short, the Meissner effect is just a natural consequence of the condensation of Cooper pairs, which are (more or less) bosons. What is important for the BCS theory is the mean-field treatment of the Cooper effect, which transmute fermions into bosons.
Please tell me if you need more details. I believe the BCS paper is really pedagogical, but please have no hesitation to ask for further details.