The BCS wave function you write down depends on a parameter $\phi$, but the ground state energy is independent of it. This implies that $\phi$ is the (would be) Goldstone mode that governs the low energy dynamics of the system. The gradient of $\phi$
is the conserved $U(1)$ current $\vec\jmath \sim \vec\nabla\phi$, and the ordinary charged current is $\vec\jmath_s =n_s e \vec\nabla\phi /m$, where $n_s$ is the superfluid density of electrons.
Because $\phi$ is a Goldstone mode the effective low energy action can only depend on gradients of $\phi$. By gauge invariance the effective action
is of the form $S[A_\mu-e\nabla_\mu\phi]$. The explicit form of $S$ can be computed from the BCS wave function, or more easily determined using diagrammatic methods. For our purposes the only important point is that $S$ has at least a local minimum if the field vanishes. This means that solutions of the classical equation of motion are of the form $A_\mu=e\nabla_\mu \phi$ (This is the London equation). Let's consider an applied electric field $\vec{E}=-\vec\nabla A_0$. I find
$$
\vec{E}=-e\vec\nabla \dot\phi=\frac{m}{n_s}\frac{d\vec\jmath}{dt}
$$
which shows that a static current corresponds to zero field, and the resisitivity is zero.
The effective action also governs other properties of the system, such as the Meissner effect, the critical current, and fluctuations of the current in a thermal ensemble.
Postscript: A commenter argues that I really need to show that $S$ has a minimum
$$
S \sim \gamma (A-e\nabla\phi)^2 + \ldots
$$
First, note that $\gamma$ determines the Meissner mass, so even without a calculation I have shown that the Meissner effect implies superconductivity. Beyond that I do indeed have to do a calculation of $\gamma$ based on the BCS wave function (I could appeal to the Landau-Ginzburg functional, but this only shifts the question to the gradient term in the LG functional). Fortunately, the calculation is straightforward, and can be found in many text books. For people with more of a particle physics interest there is a beautiful explanation in Vol II of Weinberg's QFT book. There is Anderson's famous paper on gauge invariance and the Higgs effect. I provided a version of the calculation in Sect. 3.4 of these lecture notes https://arxiv.org/abs/nucl-th/0609075
Post-Postscript: How is this different from a weakly interacting electron gas? In the electron gas I have a low energy description in terms of electrons and phonons (and other degrees of freedom). For simplicity consider the high temperature limit, where a classical description applies (as explained by Landau Fermi liquid theory, this generalizes to low T). The equation of motion for a single electron is just $m\dot v=e E$, which superficially looks like the London equation. However, this is not a macroscopic current. When I pass from microscopic to macroscopic equations there is no symmetry that forbids the appearance of dissipative terms, so the conductivity is non-zero. There is indeed a subtlety in the coupling of electrons and phonons, because without the umklapp process momentum conservation would force the conductivity to vanish.
In a superconductor the gradient of the Goldstone boson automatically describes a macroscopic current ($\gamma$ is proportional to the density of electrons). S is a quantum effective action, and dissipative terms are automatically forbidden. At finite temperature things do get a little more complicated because the total current is in general the sum of a non-dissipative supercurrent, governed by $S$, and a dissipative normal current. However, below $T_c$ part of the response is carried by a supercurrent.
