Bit of a rant coming up (and some notes that I'd like to preserve)...

If you have a process based model (i.e. one that embodies knowledge of physical, chemical and biological processes in a deterministic way, rather than stochastic way), that has many interacting, non-linear processes, how do you test for numerical stability? This kind of model is capable of exhibiting chaotic behaviour, so a small change in initial conditions may result in a large change in final output (to be truly chaotic this change need to be exponentially fast). So how do you test numerical stability? You can't change an input by a small amount and get the same or similar answer. You can't fiddle a few numerical parameters and get the same answer, as you may have changed the initial (or some other) condition, so any differences may be actual, physical, differences that are manifesting due to the new parameters including extra processes (i.e. increasing resolution would increase detail and hence what on a low resolution run is a straight line, might be a gyre on high res).

I've tried to run the same model on different architectures (windows and solaris), but I get different answers, even in double precision. But after looking deeper, this is expected. Even without any optimisation, windows (actually, intel processors) do FP calculations in 80-bit precision (as they have an 80-bit register), where possible. Sun processors do not. My thinking was that if I get the same answer on different machines, then my confidence that the model is numerically OK increases. No such luck.

However, running the code with different compilers on the same processor may be more fruitful. I'll have to see.

I must look up how GCM's are validated; they must have similar issues due to chaotic weather patterns, etc. There's a fun way to spend an afernoon!