A GISS ModelE code fragment
The several Global Climate Model or General Circulation Model (GCM) codes around the world are the work horses for the ‘what-if’ studies relative to the Global Average Temperature (and many, many other things) that appear in the IPCC Assessment Reports. There are about 20 of these codes, more or less. The actual models in the codes vary considerably from code-to-code and as a function of time with a given code. The GISS ModelE code is the NASA versions of GCMs.
Some GISS/NASA information about the ModelE GCM code is available here. There you will find information about the documentation of the model and code along with access to download the source and and an online HTML-based code viewer.
From time-to-time I like to browser around in the coding looking at properties such as in-line documentation, physical models that I can relate to, and searching for things like a data dictionary. Basically there is a notable lack of detailed documentation for either the theoretical models, numerical solution methods, and especially the contents and purpose and structure of the routines in the code. The data dictionary would be very helpful relative to understanding the variable names and units, along with a wealth of additional information. A lot of reverse engineering is required in order to understand all but the very simplest routines and functions. And, one can not be sure that the correct understanding has been extracted by use of reverse engineering.
The last time that I was in the coding I ran across the following code in subroutine DYNAM listed here:
C**** This fix adjusts thermal energy to conserve total energy TE=KE+PE
C**** Currently energy is put in uniformly weighted by mass
ediff=(TE-TE0)/((PSF-PMTOP)*SHA*mb2kg) ! C
!$OMP PARALLEL DO PRIVATE (L)
T(:,:,L)=T(:,:,L)-ediff/PK( L,:,: )
!$OMP END PARALLEL DO
starting at line 234 in the listing that I copied from the code viewer. The introductory comment is enlightening in itself. I think that this coding cannot be associated with any first-principle statement like the conservation of energy. If that was true I think the comment would have made that very clear; instead it is characterized as a ‘fix’.
The variable ‘ediff’ is the difference between the most recently-calculated sum of kinetic plus potential energy and the initial value of the sum calculated at the start of the calculation.
It has been difficult for me to track down the units of all the variables in the coding, and I have not been successful. I have determined that the temperature is very likely the potential temperature and that the variable ‘PEJ(:)’ has units [PEJ(:)] = [J/m2]; as stated in subroutine conserv_PE. I find the units situation to be kind of messy. I don’t know why the code internally does not use a purer form of SI and handle all the conversions at input and output. I suspect that eventually an error in the units will be discovered. I think the ‘! C’ at the end of the line where ‘ediff’ is calculated might mean degrees C, but I am not certain, and I think that is not correct. The coding within the ‘!$OMP PARALLEL DO PRIVATE (L)’ loop in which ediff is subtracted from ‘T(:,:,L)’ divides ‘ediff’ by another pressure, it looks like. And that pressure does not seem to be the same as in the line where ‘ediff’ is calculated.
What this ‘fix’ does, in my opinion, is adjust the calculated temperature for the difference between the recently calculated sum of kinetic plus potential energy and the initial value of the sum when the calculation was started. This might be argued to be some kind of dissipation, but that is not correct. The important point is that this adjustment is purely an artifact of the numerical methods used in the codes. The quantity ‘ediff’ is solely a result of the numerical methods and can not be based on any kind of physical-basis argument.
Thus the potential temperature reported by the NASA GISS code seems to be the potential temperature plus an artificial numerical contribution.
I do not find any diagnostic evaluation of the magnitude of ‘ediff’. If the temperature is being adjusted by an artifact of the numerical methods, I think it should be monitored. Most importantly, it should never be done in the first place.
All corrections and additional information are appreciated.