To what extent can the superconducting order parameter be thought of as a macroscopic wavefunction? I know that the order parameter does not obey the Schrodinger equation; it instead obeys the Ginzburg-Landau equation. However, I am unclear as to the situations under which the view of the superconducting order parameter as a macroscopic wavefunction breaks down, as it works to solve many problems. Why does it work so well in so many cases (e.g. in explaining flux quantization, London's equation, Meissner Effect, Josephson Effect etc.)? Where does it fail?
 A: To answer your question, one needs to understand a bit what is the Ginzburg-Landau (GL) formalism. Let us first recall the GL functional: 
$$F=\int dV\left[g\left|\left(\nabla-\dfrac{2\mathbf{i}e}{\hbar}A\right)\Psi\right|^{2}+a\left(T-T_{c}\right)\left|\Psi\right|^{2}+b\left|\Psi\right|^{4}+\dfrac{\left(\nabla\times A\right)^{2}}{2\mu}\right]$$
with $\Psi$ the (complex) order parameter, $a$, $b$ and $g$ some parameters, $\mu$ the magnetic permeability of the compounds, $V$ its volume and $T$ the temperature. Next, $T_{c}$ is the critical temperature below which the coefficient $a\left(T-T_{c}\right)$ is negative, giving rise to a finite equilibrium 
$$\dfrac{\delta F}{\delta\Psi^{\ast}}=0\Rightarrow\left|\Psi\right|^{2}=\dfrac{-2b}{a\left(T-T_{c}\right)}$$
order parameter, whereas $\Psi=0$ for $T>T_{c}$. So a short answer to your question is : whenever the GL functional is a good approximation for a superconductor will the description of this superconductor follow the GL formalism. Moreover, the GL functional is a correct description of the superconductor in the so-called GL-regime ! It sounds really like a tautology, but think about that one more time. Where are the coefficients $a\left(T-T_{c}\right)$ and $b$ coming from ? Why are there no higher-order terms, like a $c\left|\Psi\right|^{6}$ in the expansion ? What is $g$ and why is the vector-potential $A$ in the gradient term, like for charged particles ? and so on ... 
Now, let me rephrase your question as: what is $\Psi$ ? It is called an order-parameter, because it corresponds to the minimisation of the GL functional ... sounds weird to you ? It's nevertheless the definition of the order-parameter in the so-called Landau-paradigm of condensed-matter: an order-parameter is an observable quantity which is zero above some critical temperature (say $T_{c}$) and non-zero below this critical temperature (temperature can be replaced by whatever thermodynamics you want: volume, pressure, entropy, ...). I just showed that $\Psi$ follows this rule, so it is an order-parameter. 
Next question: is $\Psi$ a wave-function ? Certainly not ! The wave-function $\Phi$ of the condensate verifies the secular equation $H\Phi=E\Phi$. Bardeen, Cooper and Schrieffer (BCS) managed to give a variational wave-function of the ground-state of the Fröhlich/BCS Hamiltonian (yes, the world is full of tautology, it's clear that BCS didn't called their Hamiltonian "BCS" ...) as 
$$\Phi=\prod_{k}\left[u_{k}+v_{k}^{\ast}c_{k}^{\dagger}c_{-k}^{\dagger}\right]\Phi_{0}$$
which sounds hard to proof using the GL formalism. In fact you cannot do that in that way. 
What you can do is to follow 

L.P. Gor'kov, Microscopic derivation of the Ginzburg-Landau equations in the theory of superconductivity, Sov. Phys. JETP 9, 1364-1367 (1959). The details of the calculation is reprinted in A. A. Abrikosov, L. P. Gor’kov, and I. E. Dzyaloshinsky, Methods of Quantum Field Theory in Statistical Physics, Prentice Hall (1963).

which in brief did the following thing: 


*

*find the spectrum of the Fröhlich/BCS Hamiltonian in the mean-field approximation (Gor'kov introduced the famous anomalous Green's functions to do so)

*treat the superconducting gap self-consistently. The superconducting gap is defined as 
$$\Psi\propto\left\langle c_{k}^{\dagger}c_{-k}^{\dagger}\right\rangle $$
(compare with the above BCS Ansatz for the wave-function, note that the average $\left\langle \cdots\right\rangle \sim\left\langle \Phi\right|\cdots\left|\Phi\right\rangle $ is an average over the non-trivial (Cooper paired) vacuum, which should look like $\Phi$).


*

*expand this relation for a small gap parameter


At the end of this laborious calculation, you end up exactly with the GL functional, and you know the microscopic origins of the coefficients there. 
So, now we are able to answer the tautologic question: what is the GL-regime in which the GL-functional describes a superconductors? Well, it is clear from the Gor'kov calculation that you need $\Psi$ to be small. From the general curve of the $\Psi\left(T\right)$, you see that $\Psi$ is small when the temperature $T$ is close to the critical one $T_{c}$. In fact, a second order phase transition exhibits a linear slope of $\Psi\left(T\right)$ for $T\approx T_{c}$, which roughly gives the validity of the GL functional. 
Remark 1) Do you need the Gor'kov argument to say where the GL functional is correct ? Of course not. Even, the GL functional appeared 10 years before the Gor'kov's calculations. All you need is a bit of intuition (the order parameter is small and thus can be expanded in series because it's a second-order phase-transition, odd-powers can not be present for symmetry reason since the compounds possess an inversion symmetry, the phase should be spatially homogeneous: the penalty $g$ must be present in the functional ...) and you end up with the GL-functional.
Remark 2) It is clear that the linearised version (when $b\rightarrow 0$) of the GL functional leads to the Schrödinger equation, and so people call $\Psi$ a wave-function in that sense. One should nevertheless prefer to call $\Psi$ a complex order parameter for a charged bosonic field.
Remark 3) The GL functional is sometimes called an Abelian Higgs model, since it demonstrates the phase-fixing due to the penalty $g$ (called the condensate stiffness) in the expansion. This is entirely due to the potential vector $A$ in the GL functional, which itself is related to the fact that $\Psi$ represents a charged-bosons order-parameter. Since this model exhibits all the electrodynamics of the superconducting phase (obviously, there is no spin-dynamics in the GL functional I gave to you, so you have to adapt it if you want to describe the spin-dynamics as well), the GL functional is described as the model to discuss superconductivity, without introducing it as a generalisation of the Landau model for phase-transition and order-parameter, here adapted to the case of a charged quantum-field. I guess the origin of this question lies in the absence of knowledge on the Landau model of second-order phase transition. Also, to introduce the GL functional without microscopic discussion is often done in high-energy textbooks which want to focus on the Higgs-mechanism.
Remark 4) Can we apply the GL functional to describe the superconductors at low temperatures ? In principle not, but there are some arguments showing that the Abelian Higgs model can be generalised to low temperatures, see e.g. Greiter: Is electromagnetic gauge invariance spontaneously violated in superconductors?
I realised after writing this long answer that I haven't answered your last question : why the GL functional describes so well the electrodynamics of superconductors ? The reason is that the Higgs mechanism is the unique modification you need to do to pass from the electromagnetism of the normal metal to the electromagnetism of the superconductor. More heuristically, since the Cooper pairs move without resistance, it is clear that they will make-the-law in disordered systems. So you only need to capture their fundamental property, which is the pure diamagnetism, see e.g. this question elsewhere on this website (note that Meissner effect can be replaced by perfect diamagnetism all the answers to this question).
