Multi-scale/level computer models are efficient for simulating solutions of partial differential equations on nontrivial geometries using a large number of modes. This is because fine-scale structures in the quantities of interest cannot, in practice, be resolved using a single-level finite element method (FEM) applied directly to the full model, or standard spectral basis functions are not appropriate for the nontrivial geometry and related constraints. The multi-scale/level FEM (MFEM) is based on first applying a FEM for a reduced model problem, obtained by restriction to either (i) several local subdomains or (ii) a dominant linear differential operator in the model. The resulting structure-aware precomputed solutions are then used as basis functions to simulate the full model in the MFEM framework. If these are eigenfunctions of a reduced model eigenvalue problem, such basis functions provide spectrally accurate approximations to the solution of the full model and the resulting method is called the spectral MFEM (SMFEM). In this article we develop, analyze, and implement time-space fully discrete implicit SMFEM (ISMFEM) algorithms for efficiently simulating a Ginzburg--Landau (GL) system modeling superconductivity on a class of superconducting surfaces $S$. The spatial linear differential operator in the time-dependent nonlinear GL model is the Schrödinger operator $( i \nabla + {A}_0)^2$ on $S$, with magnetic vector potential ${A}_0$. For our ISMFEM, the spectrally accurate basis functions are eigenfunctions of a reduced stationary linear model governed by $(i \nabla + {A}_0)^2$. In a recent work we developed, analyzed, and implemented an algorithm to compute these functions using an analytic-numeric hybrid method, in conjunction with a high-order FEM. Using these precomputed solutions as basis functions, the spectrally accurate algorithms developed and analyzed for our new GL computer models on nontrivial geometries are validated using both nonlinear and linearized implicit time discretizations.