I finally have good versions of the code, which have passed every test I know!
One version is very slow, but very accurate. It runs with a time step of 0.9 seconds only, which gives us a runtime of 325 days for a simulation of 11 years i.e. one sunspot cycle. This code can be used for observing small scale features which do not exist for long.
The faster version of the code runs with a time step of 250 seconds, which gives us a runtime of 24 hours per sunspot cycle. In this code, the grid points within 4.32 degrees of latitude from the pole have 50 points in the toroidal direction, whereas the rest of the latitudes have 1000 points in the toroidal direction.
The evolution of the radial field is second order accurate in both forms of the code, and there is a slight flux imbalance due to numerical round off, which is less than 0.1 % of the total flux.
I ran a simulation with the modeled sunspot input for one solar cycle, and the output is shown in the video here.