Monday, November 30, 2009

Was the CRU putting their Thumb on the Scale?

I downloaded the CRU data and glanced at a couple of files. One that I looked at was: documents/cru-code/f77/mnew/sh2sp_normal.for. A comment at the top explained that this subprogram "converts sun hours monthly .nrm to sunshine percent (n/N)". I glanced at the code and noticed this:

DATA DYS/31,28.25,31,30,31,30,31,31,30,31,30,31/

They're going to take the number of hours of sunshine in a month and divide it by the number of days in that month.

That quarter day in February didn't sit well with me. I'm going to speculate that in the real world where hours of sunshine are observed, they are observed (or not) in whole 24 hour days. In other words, you can't use the coding "trick" of slathering the quadrennial leap day over every year—ya gotta put it in the year it's observed. I do sympathize. Accurate calendar accounting is much more difficult than is generally realized.

The good news is that it doesn't matter. Let's look at a little more code, so I can explain why:

integer dys(12)
real day(12),dec(12),sunp(12)
DATA DYS/31,28.25,31,30,31,30,31,31,30,31,30,31/

We can think of that 28.25 days as being a square peg and that integer array, dys, as being a round hole. When we assign a real number to an integer variable, Fortran helpfully truncates that pesky decimal portion, leaving just 28 in our case. That means our leap day doesn't even show up every four years.

Fortran is not case sensitive, so both DYS and dys are the same variable. I think switching case as they've done above is poor form. It's certainly confusing to me, but it's been awhile since I coded in Fortran.

Luckily, there's something you can do to help. You can buy Bad Code Offsets for the East Anglia Climate Research Unit to mitigate the harm done by their bad code! (H/T: L'Ombre de o'livier)

No comments: