$3 \times 3 $ Magic Square of Squares

You will want to look at A search for $3\times3$ magic squares having more than six square integers among their nine distinct integers, by Christian Boyer, and at the papers by Bremner and others that Boyer references. You can't hope to find a new one until you understand the methods used to find the one that's already known.

Boyer also published a paper which I haven't seen: Some notes on the magic squares of squares problem, Math. Intelligencer 27 (2005), no. 2, 52–64, MR2156534 (2006d:05024).

I've dealt with this some monthes ago and have this in some scribbles. Don't know whether this is of any help.
From the ansatz (where m is the horizontal,vertical and diagonal sum) $$ \begin{array} {rrr|r} & & & m \\ a^2 & b^2 & c^2 & m \\ d^2 & e^2 & f^2 & m \\ g^2 & h^2 & i^2 & m \\ \hline m&m&m&m \end{array} $$ writing this as an array of equations and using Gauss-reduction I arrived at the following magic square with three free parameters e,i,h $$ \begin{array} {rrr|r} & & & 3e^2 \\ 2e^2-i^2 & 2e^2 -h^2& -e^2 +h^2 + i^2 & 3e^2 \\ -2e^2+h^2+2i^2 & e^2 & 4e^2-h^2-2i^2 & 3e^2 \\ 3e^2-h^2-i^2 & h^2 & i^2 & 3e^2 \\ \hline 3e^2&3e^2&3e^2&3e^2 \end{array} $$ and if I recall correctly I've seen, that the parameters must be odd and not divisible by 3 or in other words $e^2,i^2,h^2$ must be congruent 1 modulo 12 .
I didn't proceed then, however perhaps that representation is of some interest for you.

[update2] Here is one more information which I forgot to include earlier. We can express the conditions on the entries in the magic square, which should also be squares, depending on the three parameters $e,h,i$ as a small matrix-multiplication: $$ \begin{array} {r} &&&&|&e^2| \\ &&&&|&h^2| \\ &&&&*|&i^2| \\ \hline |& 2 & 0 & -1| & |&a^2| \\ |& 2 & -1 & 0| & |&b^2| \\ |&-1 & 1 & 1| & = |&c^2| \\ |&-2 & 1 & 2| & |&d^2| \\ |& 4 & -1 &-2| & |&f^2| \\ |& 3 & -1 &-1| & |&g^2| \\ \end{array} $$ I found that mucht tempting to try to make something out of that structural description to say something about the possibilities for all entries simultanously to be squares, but have not yet a better expression.

[update] The comments below motivated me to simply try out that parametrized problem. Using a Pari/GP-routine with a three-fold loop for the base-parameters $e,h,i$ I got this 7-square-solutions in 53 secs: (which is also the given 7-squares solution shown in the thread's initial question) $$ \small \begin{bmatrix} 205^2 & 527^2 & 222121 \\ 360721 & 425^2 & 23^2 \\ 373^2 & 289^2& 565^2 \end{bmatrix} \small \begin{bmatrix} 222121 & 527^2 & 205^2 \\ 23^2 & 425^2 & 360721 \\ 565^2 & 289^2 & 373^2 \end{bmatrix} $$ The symmetry of the two solutions indicate, that I could have halved the consumption of time If I had some smarter search-criteria.

With some improved criteria for the loop (100 sec, $e$ used up to 3000) I found some more - unfortunately the're only the trivial multiples of the first solution... : $$ \small \begin{matrix} & a^2 & b^2 & c^2 & d^2 & e^2 & f^2 & g^2 & h2 & i^2 \\ \hline & 410^2 & 1054^2 & 2^2\cdot 151 \cdot 1471 & 2^2\cdot 137 \cdot 2633 & 850^2 & 46^2 & 746^2 & 578^2 & 1130^2 \\ & 615^2 & 1581^2 & 3^2\cdot 151 \cdot 1471 & 3^2\cdot 137 \cdot 2633 & 1275^2 & 69^2 & 1119^2 & 867^2 & 1695^2 \\ & 820^2 & 2108^2 & 4^2\cdot 151 \cdot 1471 & 4^2\cdot 137 \cdot 2633 & 1700^2 & 92^2 & 1492^2 & 1156^2 & 2260^2 \\ & 1025^2 & 2635^2 & 5^2\cdot 151 \cdot 1471 & 5^2\cdot 137 \cdot 2633 & 2125^2 & 115^2 & 1865^2 & 1445^2 & 2825^2 \\ & 1230^2 & 3162^2 & 6^2\cdot 151 \cdot 1471 & 6^2\cdot 137 \cdot 2633 & 2550^2 & 138^2 & 2238^2 & 1734^2 & 3390^2 \\ & \vdots \\ k^2*\ldots&205^2 & 527^2 & 151 \cdot 1471 & 137 \cdot 2633 & 425^2& 23^2 & 373^2 & 289^2 & 565^2\\ & \vdots \end{matrix} $$ Obviously there is no number $k^2$ which would make the entries in columns $c^2$ and $d^2$ a perfect square, so this scheme cannot provide a better solution for higher $k$.

Here is the Pari/GP-code (updated)

isin(x,vgl)=if (x<1,return(1)); for(k=1,#vgl,if(x==vgl[k],return(1)));return(0);

{ listsqsq(max_e=100,max_nosq=3,min_e=1)= local(a,b,c,d ,f,g,  no_sq,a2,b2,c2,d2,e2,f2,g2,h2,i2,list,li);
 list=vectorv(20000);li=0;
 for(e=min_e,max_e, e2=e^2;
  for(h=1,ceil(1.5*e),  if(h==e,next()); \\ no higher h needed
        h2=h^2;
        b2=2*e2 - h2; if(isin(b2,[e2,h2]), next()); 
             if(issquare(b2)==0, next()); 
             b=sqrtint(b2);

        for(i=sqrtint(e2-ceil(h2/2)),sqrtint(2*e2-floor(h2/2))+1,  if(isin(i,[e,h,b]),next());  \\ no higher i needed
              i2=i^2;no_sq=0;
              g2=   b2 +e2-i2;  if (isin(g2,[b2,e2,h2,i2])            ,next()); if(issquare(g2)==0,next()); g=sqrtint(g2);
              a2= 2*e2 -i2;     if (isin(a2,[b2,e2,h2,i2,g2])         ,next()); if(issquare(a2)==0,a=-a2;no_sq++, a=sqrtint(a2));
              c2=-e2+h2+i2;     if (isin(c2,[b2,e2,h2,i2,g2,a2])      ,next()); if(issquare(c2)==0,c=-c2;no_sq++, c=sqrtint(c2));
              d2=-2*e2+h2+2*i2; if (isin(d2,[b2,e2,h2,i2,g2,a2,c2])   ,next()); if(issquare(d2)==0,d=-d2;no_sq++, d=sqrtint(d2));
              f2= 4*e2-h2-2*i2; if (isin(f2,[b2,e2,h2,i2,g2,a2,c2,d2]),next()); if(issquare(f2)==0,f=-f2;no_sq++, f=sqrtint(f2));
              if(no_sq>max_nosq,next());
              idx=prime(e)*prime(h)*prime(i)*prime(g);
              li++;list[li]=[a,b,c,d,e,f,g,h,i,log(idx),no_sq];
            )
      );
    );
   if(li==0  , return(Mat([0])));
   list=Mat(vecextract(list,Str("1..",li)));
   list=vecsort(list~,[10,4,7,8,3,6,7])~;
   return(list);}

This isn't a solution. But it's too big to be a comment. I just hope it's useful.

First observation

You can always arrange the numbers in a magic square so that the smallest is top-middle and the next is bottom-right, then there are two possible arrays of ranks of the $9$ numbers.

 TYPE 1        TYPE 0
 8  1  6       8  1  7
 3  5  7  and  4  5  6
 4  9  2       3  9  2

The first is itself a magic square, the second is not. The nicest example of a TYPE $0$ magic square is

 8   0  7
 4   5  6
 3  10  2

Second Observation

For all magic squares, the sum of the squares of the first row (column) equals the sum of the squares of the last row (column). For example

 8² + 0² + 7² = 3² + 10² + 2² = 113

So you are not only looking for numbers such that $$a+b+c=A+B+C$$

and

$$a^2 + b^2 + c^2 = A^2 + B^2 + C^2$$

you must also have

$$a^4 + b^4 + c^4 = A^4 + B^4 + C^4$$

Even further, $\{a,b,c\}$ and $\{A,B,C\}$ must both be complete residue systems modulo $3$.