The Day The Thermometer Music Died. Thermometers by Year Crashes.
And they were singin’: “Bye Bye Miss American Pie, Drove my GIStemp to the Levy But the Levy was Dry; and them Good Ol’ Boys was Drinking Whiskey and Rye, and Singin’ ‘This Will Be The Day That I Die!, This Will Be The Day, That I Die’…”
Introduction to GHCN – The Global Historic Climate Network
This is an ‘aggregator posting’. I’m putting here the links to the various individual analysis steps of GHCN on a global basis. (This is so that, in the future, I only need to use this one link to get to any of them).
So what are the postings?
I looked at GHCN input data from various places around the world. By continent. By major country on some continents. As time permits, I’ll add more fine grained looks at some other countries. (Under the Asia thread, I found that Japan now has no thermometer above 300 M. Who knew Japan was as flat as Kansas… So I’m going to “do Japan” at some time and see what else turns up… When that happens, a link will be added here.) With that, here is the list of links to “What the GHCN (Global Historic Climate Network) data look like, by continent and with selected countries”.
GHCN is the Global record of land thermometers (that is the “historic” part – clearly their bias is that satellites are the future and those actual instruments on the ground are so ‘historic’ as to be positively ‘old school’). Frankly, a well tended mercury thermometer is hard to beat, but I’m not in charge (and they are not well tended
http://www.surfacestations.org ).
You can get a bit more detail, along with some file format information at:
http://chiefio.wordpress.com/2009/02/24/ghcn-global-historical-climate-network/
You can download your own copies of the GHCN data from the ftp site at:
ftp://ftp.ncdc.noaa.gov/pub/data/ghcn/v2
Which also includes more information in “readme” files.
The Links to my Analysis
These are presented in an order different from their original writing. You may well find that when you read one it talks about “what we saw before” that you have not yet seen. “No worries”. There is not a lot of time dependence and you ought to be able to take them in any order without too much of an ‘issue’.
The World, and The Hemispheres
Early on, I noticed that the history of the thermometer record had “The Thermometers March South”. Initially I assumed this was just an artifact of the spread of technology, and time, spreading from the north to the south. And perhaps some spread of wealth and thermometers in the Jet Age as airports spread to tropical vacation lands. Yet there was an odd discontinuity at the end. In the 1990’s, the thermometer count plunged overall, and the percentage in the Northern Cold band was cut dramatically. This was the early investigation that lead to all the other links here:
http://chiefio.wordpress.com/2009/08/17/thermometer-years-by-latitude-warm-globe/
Early on, before I had worked out the details and polished the code enough to take detailed slices through the GHCN data; I had made some postings that looked at large swathes of the data. The selections were often done by hand (using Linux / Unix tools like ‘grep’) and the processes later were turned into the scripts that are listed at the bottom of this posting.
One early cut looked at the Northern Hemisphere in total. In retrospect, you can see the “Day the Thermometer Music Died” in the charts of winter months. Thermometers at the beach in Los Angeles don’t get very cold in winter, those at Squaw Valley Ski Area do. In California, all our thermometers have left the mountains and are now on the beach with 3/4 of them near L.A. and San Diego. I didn’t know that when this posting was made, but you can clearly see the effect of these changes in the Northern Hemisphere graphs:
http://chiefio.wordpress.com/2009/10/18/the-northern-hemisphere-what-warming/
I initially thought that the stability in the record seen in long lived thermometers must have been an artifact of modification flag changes with the changes of equipment (as many placed moved to automated temperature recording). Since then, I’ve determined that there really were a bunch of thermometers deleted. 90+% in many major countries around the world. There is still a discontinuity at the end, and that lead to more detailed investigations.
Sidebar on Data Sources: In some of the links, you will see the file v2.mean as the data source. That is the straight GHCN data. In others, you will see v2.mean_comb as the data source. This is the result of “STEP0″ of GIStemp. It has had the Antarctic data enhanced with some extra data from the Antarctic projects directly; it has had an extended copy of one site in Germany used to replace the original; and it has had the USHCN copy of the US data merged with the GHCN copy (and that will only effect the US data).
Because of this, for substantially everywhere in the world other than the USA and Antarctica, these two data sets will have indistinguishable effects.
For Antarctica it is valuable to use the v2.mean_comb file to see most of the place (though a look with only the v2.mean might be instructive about GHCN…); while for the USA, since GIStemp did not make the transition to USHCN.v2 (version 2) the impact of the USHCN data ‘cuts off’ in 2007. There is no difference between the GHCN data for the USA and the v2.mean_comb contents after that date.
The Detail Study Links
We will step through the whole globe, region by region.
North America:
This first one has California in the title, and there is an important issue about California in particular in the posting; but the posting is about all of The U.S.A. and uses data for the whole country. (The “issue” is that as of 2009, California has 4 active thermometers in GHCN. 3 on on the beach in Southern California and one is at the airport in San Francisco, we presume waiting for it’s ride to L.A….)
http://chiefio.wordpress.com/2009/10/24/ghcn-california-on-the-beach-who-needs-snow/
This posting looks at changes in Mexico (which is found to have a strong thermometer change bias) but also looks at the “little bits” left over in North America when you leave out Canada, the USA, and Mexico. Thermometers can’t move very far in Belize or The Bahamas… and we find that the temperature record there shows no warming (and even small hints of cooling).
http://chiefio.wordpress.com/2009/11/01/ghcn-mexico-a-megathermal-vacation-band/
The Arctic:
In this posting we look at the Canadian and Russian arc that surrounds most of the Arctic. Russia is split into a “European” and “Asian” chunk in GHCN anyway, so I hope this geographic discontinuity is not too jarring… But I think it does make sense to cover Canada under an “arctic” listing.
http://chiefio.wordpress.com/2009/10/27/ghcn-up-north-blame-canada-comrade/
Similarly, we put the Nordic Europe in one bucket to see what it looks like in isolation.
http://chiefio.wordpress.com/2009/10/29/nordic-north-nothing-much-to-see/
Once again we see that when the thermometer record in a country is grossly distorted by deletions, there is an artifact (though not always in the direction expected) and that when a geography is instrumented with a more stable thermometer set, there is no warming present.
Basically, if we’re setting up a global calorimeter to measure heat gain / loss, we need to stop changing the instrument by moving around thermometers and adding / deleting them. Pulling 90% of your thermometers out of the calorimeter makes calibration impossible. (And that renders the results more fantasy than useful. Heck, it makes “Cold Fusion” calorimetry look positively stellar in comparison…)
The Pacific Islands and Australia / New Zealand:
http://chiefio.wordpress.com/2009/10/29/ghcn-pacific-basin-lies-statistics-and-australia/
http://chiefio.wordpress.com/2009/10/23/gistemp-aussy-fair-go-and-far-gone/
And one of my favorites where we see how one island can shift the whole region:
http://chiefio.wordpress.com/2009/11/01/new-zealand-polynesian-polarphobia/
The end of it all is that the entire Pacific Basin is substantially flat on temperatures. Hard to have “Global Warming” if the Pacific is not participating. Australia and New Zealand show warming, but only due to thermometer change artifacts. For New Zealand, it is one single cold thermometer: And when that one is deleted from the whole record, not just the last few years, New Zealand has no “Global Warming” either.
Hard to have “Global Warming” when the 1/2 of the planet that is the Pacific Basin is dead flat with only a small “ripple” as the PDO flips state every 30 or so years.
The Antarctic (we covered one pole, let’s do the other):
One of the more interesting bits is in an update way down at the bottom. I broke out bits of Antarctica by geography so folks can compare east to west and peninsula to center. One site shows 4 years of dirty data, but the NASA site GIStemp map somehow turns this into one single data point; though in the wrong year! Another site has dead flat 21.x C entry and exit to the data series, yet the NASA GIStemp chart has the entry and exit ends of the graph flip flopping like a fish on the dock by about 2 C. (Yes, 2 whole degrees, forget the tenths place…) This, IMHO, is clear proof that the GIStemp process and NASA charts are as much fantasy as anything else.
http://chiefio.wordpress.com/2009/11/02/ghcn-antarctica-ice-on-the-rocks/
South America:
Hard to have “Global Warming” when most of South America is not participating…
http://chiefio.wordpress.com/2009/11/02/ghcn-south-america-siesta/
http://chiefio.wordpress.com/2009/10/24/ghcn-brazil-sambas-north/
http://chiefio.wordpress.com/2009/10/24/argentina-cool-on-the-pampas/
Africa:
Where we find that the continent is not warming, though the thermometer coverage moves around quite a bit, we can still see everything we need to see:
http://chiefio.wordpress.com/2009/10/29/africa-halle-barry-hot-steady-with-variable-coverage/
Africa is a hot place, but it is not getting hotter. Hard to have “Global Warming” when Africa is not participating…
Asia:
The bulk of the countries in Asia have no “Global Warming”. They are a fairly smooth temperature set. It is only the Siberian thermometer changes and the Chinese thermometer deletions that show much change. There are hints in the “Without Siberia and China” charts of other things to dig into, but the “Global Warming” signal is definitely squashed by taking out the thermometer “issues” in the two big countries. This “hint”, though, led to an early look at altitude changes in Japan, where we find it no longer has any thermometers over 300 meters elevation. It seems that, like California, Japanese thermometers like it on the beach…
http://chiefio.wordpress.com/2009/11/02/ghcn-asia-chinese-footprints-in-siberian-snow/
China, too, has had a thermometer count crash:
http://chiefio.wordpress.com/2009/10/28/ghcn-china-the-dragon-ate-my-thermometers/
Europe:
http://chiefio.wordpress.com/2009/11/02/ghcn-europe-goes-mediterranean/
Europe is an interesting study in that it is one of the earlier places where thermometer migration south shows up. It has a smoother character to the change. It is harder to see the time and spacial onset, and harder to pick a “smoking gun” moment; but it is a strong example of The March Of The Thermometers southward.
GEEK Corner: The Computer Code
I will also be putting here the listings of the code I used to process the GHCN data. This is so that anyone who wishes to duplicate any of this work can see what I did ( and hopefully both replicate it and improve on it).
This code is written in “bash” for the scripts (that ought to run in ’sh’ and ‘ksh’ environments too, if I did things right) and in FORTRAN for the main code. Why FORTRAN? Simply because I’m ‘deconstructing GIStemp’ and it is written in FORTRAN. I find it easier to only have one language loaded into my brain at a time (well, 2 if you count bash… 3 if you include English… but I use Spanish at the fast food places… and French is currently on line too… and… well, lets just say it’s crowded in here and the less added stuff the better). Besides, C is somewhat antithetical to FORTRAN file formats and these data come from FORTRAN programs. So it’s easier to just stick with “the horse what brung you to the party”…
If you don’t like it, perhaps a plea to the Gods Of Source Code will result in someone posting a C translation. Other than the file format issues, it is a trivial bit of code to produce.
One Disclaimer: All this code was written as a fast “hand tool” cobbled together from other code as a base. It isn’t the best solution, only the more expedient one. There is plenty of room for improvement…
Finally, the scripts may extend off the right side of the page. WordPress, in this theme, gives me the Hobson’s Choice of doing a “preformatted” listing and truncating visibility on the right, or not, and letting it steal all the white space and wipe out the formatting. For those programmers who really want to see the “stuff off the right edge” you can just choose “view page source” for your browser and all the text is there. Programmer types ought not to have much of a problem with that, and non-programmers won’t care that they don’t see the right most text.
Merge Temperature History with Station Information
First, we make a merged file from the v2.mean and v2.inv files. Horridly inefficient, but I was more interested in getting done quickly than elegance. It ought to be a straight database load, but you do what is fastest to complete some times.
[chiefio@tubularbells vetted]$ cat ccaddlat.f
C Program: ccaddlat.f
C Written: Nov 3, 2009
C Author: E. M. Smith
C Function: This program matches v2.inv and v2.mean on Station IDs
C sorted in order by cc stationID(8) and produces a merged v2.mean
C format file.
C So as to match station location info with thermometer records.
C
C Copyright (c) 2009
C This program is free software; you can redistribute it and/or modify
C it under the terms of the GNU General Public License as published by
C the Free Software Foundation; either version 2, or (at your option)
C any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C General Housekeeping. Declare and initialize variables.
C
C There are some left in there from a prior version of the code.
C I've not "preened" this version for consistency of definitins
C declarations and use. -ems
C
C itmp - an array of 12 monthly average temperatures for any
C given record. The GHCN average temp for that station/year.
C icc - the "country code"
C iyr - The year for a given set of monthly data (temp averages).
C idcount - the count of valid data items falling in a given month
C for all years. An array of months-counts with valid data.
C sid - Station ID for the idcount.
C line - A text buffer to hold the file name of the input file
C to be processed, passed in as an aguement at run time.
C line2 - A text buffer for the file name of Station IDs to process.
C oline - A text buffer for the output file name
C cid - Country Code - Station ID 3+8 char
C csid - Country Code - Station ID 3+8 char
C id - Staion ID in the input data.
C buff - A text buffer for holding input / output record data.
C buff2 - A text buffer for holding input / output record data.
C2345*7891 2 3 4 5 6 712sssssss8
integer itmp(12), icc, iyr, idcount, lc
character*128 line, line2, oline
character*8 id, sid
character*11 cid, csid
character mod
character*95 buff
character*64 buff2
data itmp /0,0,0,0,0,0,0,0,0,0,0,0/
data buff /" "/
mod=""
icc=0
id=""
sid=""
cid=""
csid=""
idcount=1
iyr=0
lc=0
C Get the name of the input file, in GHCN format. The file must be
C sorted by ID (since we count them in order.)
C The name of the output file will be inputfile.withmean
C line2 will hold the v2.mean type file sorted on ID (9 char)
C line will hold the v2.inv info sorted on ID (8 char)
call getarg(1,line)
call getarg(2,line2)
oline="v2.mean.inv"
C2345*7891 2 3 4 5 6 712sssssss8
C lines of the form "write(*,*)" are for diagnostic purposes, if desired.
C write(*,*) oline
open(1,file=line,form='formatted')
open(2,file=line2,form='formatted')
open(10,file=oline,form='formatted') ! output
C write(*,*) line2
C read in a v2.mean record, then a v2.inv description record
read(2,'(a11,a1,a64)',end=500) csid, mod, buff2
C write(*,*) "Before 20: ", icc, sid, mod, buff2
20 read(1,'(a11,a95)',end=400) cid,buff
C write(*,*) "After 20: ", icc," ", id, " ", buff
C write(*,*) "at 20: ", sid," ",mod," ",buff2
40 continue
if(cid.eq.csid) goto 70
if(cid.gt.csid) goto 80
if(cid.lt.csid) goto 90
C write(*,*) "Compiler error!. You can't get here!"
70 DO WHILE (cid .eq. csid)
write(10,'(a11,a1,a64,1x,a95)') csid,mod,buff2,buff
read(2,'(a11,a1,a64)',end=300) csid, mod, buff2
C write(*,*) "From 70: ", icc," ", sid," ", mod," ", buff2
END DO
goto 40
80 DO WHILE (cid .gt. csid)
read(2,'(a11,a1,a64)',end=300) csid, mod, buff2
C write(*,*) "From 80: ", csid," ", mod," ", buff2
lc=lc+1
if (lc.gt.1) then
write(*,*) "in 80: ", csid," ",mod," ", buff2
write(*,*) "in 80: ", cid," ", buff
write(*,*) "Too many times in loop 80", lc
end if
END DO
lc=0
goto 40
90 DO WHILE (cid .lt. csid)
read(1,'(a11,a95)',end=200) cid,buff
C write(*,*) "From 90: ", icc, id, buff
C write(*,*) "id vs sid", id, sid
END DO
goto 40
200 continue
STOP "Out of v2.inv Station Records - ought to be rare!"
300 continue
STOP "Out of v2.mean records. Most likely case."
400 STOP "Input file 1 blank on first record!"
500 STOP "Input file 2 blank on first record!"
END
[chiefio@tubularbells vetted]$
I used g95 as the FORTRAN compiler. The “wrapper script” that does the environmental set up and manages the program is somewhat overly complex. It lets you choose various input files other than the GHCN v2.mean and v2.inv files so that you can use the same tools on other sources of data. You don’t really need any of that (but it is helpful in GIStemp analysis, so I can choose to use the v2.mean.z file after Antarctica is added, or not…) Mostly it just calls the one program and hands it the v2.inv and v2.mean files as input.
[chiefio@tubularbells vetted]$ cat mkinvmean
echo " "
echo "Optional: sort the v2.inv type file by Country Code / Station ID"
echo " "
ls -l ${1-/gnuit/GIStemp/STEP0/input_files/v2.inv} v2.inv.ccid
echo " "
echo -n "Do the sort of v2.inv data into v2.inv.ccid (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
INVFILE=v2.inv.ccid
echo INVFILE= $INVFILE
sort -n -k1.1,1.11 ${1-"/gnuit/GIStemp/STEP0/input_files/v2.inv"} > $INVFILE
else
INVFILE=${1-"/gnuit/GIStemp/STEP0/input_files/v2.inv"}
echo INVFILE= $INVFILE
fi
echo " "
echo "v2.inv data from:"
echo " "
ls -l $INVFILE
echo " "
ls -l ${2-"/gnuit/GIStemp/STEP0/input_files/v2.mean"} v2.sort.ccid
echo " "
echo "Optional: Sort the v2.mean type file by Country code / station ID"
echo " "
echo -n "Do the sort of v2.mean into v2.sort.ccid (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
MEANFILE=v2.sort.ccid
sort -n -k1.1,1.11 ${2-"/gnuit/GIStemp/STEP0/input_files/v2.mean"} > $MEANFILE
else
MEANFILE=${2-"/gnuit/GIStemp/STEP0/input_files/v2.mean"}
fi
echo " "
echo "v2.mean data from: "
echo " "
ls -l $MEANFILE
echo " "
echo "Then feed the v2.inv data, sorted by cc/station ID, and "
echo "v2.mean by ccid to the program that matches ID to description"
echo "records in the data set."
echo " "
echo -n "Do the matching of $INVFILE with $MEANFILE into: ./v2.mean.inv (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
bin/ccaddlat $INVFILE $MEANFILE
fi
echo " "
echo "Produced the list of v2.mean.inv records "
echo "(temps, with station data, sorted by CCStationID"
echo " "
ls -l v2.mean.inv
echo " "
ls -l v2.sort.ccid v2.inv.ccid
echo " "
echo -n "Remove the work files v2.sort.ccid and v2.inv.ccid (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
echo rm v2.sort.ccid v2.inv.ccid
rm v2.sort.ccid v2.inv.ccid
fi
echo " "
echo -n "Does v2.mean.inv need sorting by CC, Station ID, and Year (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
# sort -n -k1.1,1.16 v2.mean.inv
echo bin/meansortidyr v2.mean.inv giving v2.mean.inv.ccsidyr
bin/meansortidyr v2.mean.inv v2.mean.inv.ccsidyr
echo
ls -l v2.mean.inv.ccsidyr
fi
[chiefio@tubularbells vetted]$
The script “bin/meansortidyr” is a one line sort that could easily be put ‘in line’ in the above script. I broke it out as a convenient ‘hand tool’:
[chiefio@tubularbells vetted]$ cat bin/meansortidyr
sort -n -k1.1,1.16 ${1-”../../STEP0/input_files/v2.mean”} > ${2-”v2.sort.ccidyr”}
[chiefio@tubularbells vetted]$
Temperature Averages By Years
This program creates a history of temperature over years. You could run it against the above combined v2.mean v2.inv file or against the v2.mean format files of GIStemp and GHCN (it only uses the first v2.mean format part of the file, so will work with any of them.)
[chiefio@tubularbells analysis]$ cat src/lmyears.f
C2345*7891 2 3 4 5 6 712sssssss8
C Program: lmyears.f
C Written: October 30, 2009
C Author: E. M. Smith
C Function: To produce a list of Global Average Temperatures for
C each year of data in a GHCN format file, with one GAT for each
C month and a total GAT for that year. Summary GAT records are
C produced for the whole data set as a "crossfoot" cross check of
C sorts. While you might think it silly to make a "global average
C temperature" for a 130 year (1880 to date) or 308 year (1701 the
C first data in GHCN, to date) interval, once you accept the idea
C of adding together 30 days, or 365 days of records, or
C thermometers from all over the planet "means something":
C Where does it end?
C Personally, I think the whole idea of a GAT is bogus,
C but if you accept it as a concept (and GIStemp and the AGW
C movement do) then you must ask:
C "in for a penny, in for a pound":
C When does the GAT cease to have some value, and exactly why?...
C
C So I produce GAT in several ranges and you can inspect it
C and ponder.
C Copyright (c) 2009
C
C This program is free software; you can redistribute it and/or
C modify it under the terms of the GNU General Public License as
C published by the Free Software Foundation; either version 2,
C or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You will notice this "ruler" periodically in the code.
C FORTRAN is position sensitive, so I use this to help me keep
C track of where the first 5 "lable numbers" can go, the 6th
C position "continuation card" character, the code positions up
C to number 72 on your "punched card" and the "card serial number"
C positions that let you sort your punched cards back into a proper
C running program if you dropped the deck. (And believe it or
C not, I used that "feature" more than once in "The Days Before
C Time And The Internet Began"...
C2345*7891 2 3 4 5 6 712sssssss8
C Oddly, within the 7-72 positions, FORTRAN is not position
C sensitive. This was so if you put in an accidental space
C character, you didn't need to repunch a whole new card...
C Oh, and if you type a line past the 72 marker, you can cut
C a variable name short, creating a new variable, that FORTAN
C will use as an implied valid variable. So having "yearc"
C run past the end can turn it into "year" "yea" "ye" "y"
C which will then be the actual variable you are using in that
C line, not the one that prints out on your card. The
C source of endless bugs and mirth 8-}
C General Housekeeping. Declare and initialize variables.
C
C itmp - an array of 12 monthly average temperatures for any
C given record. The GHCN average temp for that
C station/year.
C incount - the count of valid data items falling in a given month
C for a year. An array of months-counts with valid data.
C nncount - the count of valid data items falling in a given month
C for all years. An array of months-counts w/ valid data.
C itmptot - Array of the running total of temperatures, by month.
C ymc - count of months for the total of a year with some
C valid data.
C ntmptot - running total of all temperatures, by month column.
C icc id iyr nyr iyrmax m iyc - countrycode, Stn ID, on year
C of data as monthly averages of MIN/MAX temps, max year
C so far, month, iyc In Year Counter: # recs in year.
C tmpavg - Array of average temperatures, by month. The data
C arrive as an INTEGER with an implied decimal point
C in itmp. This is carried through to the point where
C we divide by 10 and make it a "REAL" or floating point
C number in this variable.
C ttmpavg - Total of temperature data by month for all years.
C tymc - count of months for the total of all data with some
C valid data.
C eqwt - Total of monthly averages of temperature data, by month.
C eqwtc - Counter of months with valid data.
C eqwtg - Grand Total of calculated monthly averages of MIN/MAX
C averages. Divided by eqwtc for Grand Avg.
C gavg - Global Average Temperature. GAT is calculated by
C summing tmpavg monthly averages that have valid data,
C then dividing by the count of them with vaid data.
C ggavg - The Grand Grand Average Temperature,
C whatever it means...
C line - A text buffer to hold the file name of the input file
C to be processed, passed in as an aguement at run time.
C oline - A text buffer for the output file name,
C set to the input_file_name.GAT
C2345*7891 2 3 4 5 6 712sssssss8
integer incount(12), nncount(12), itmptot(12), ntmptot(12)
integer itmp(12)
integer icc, id, iyr, nyr, iyrmax, m, iyc
real tmpavg(12), ttmpavg(12), eqwt(12), eqwtc(12)
real gavg, ggavg, ymc, tymc, eqwtg
character*128 line, oline
data incount /0,0,0,0,0,0,0,0,0,0,0,0/
data nncount /0,0,0,0,0,0,0,0,0,0,0,0/
data itmptot /0,0,0,0,0,0,0,0,0,0,0,0/
data ntmptot /0,0,0,0,0,0,0,0,0,0,0,0/
data itmp /0,0,0,0,0,0,0,0,0,0,0,0/
C data tmpavg /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./
C2345*7891 2 3 4 5 6 712sssssss8
data tmpavg /-99.,-99.,-99.,-99.,-99.,-99.,-99.,-99.
*,-99.,-99.,-99.,-99./
data ttmpavg /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./
data eqwt /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./
data eqwtc /0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0./
icc =0
id =0
iyr =0
nyr =0
iyrmax=0
m =0
iyc =0
gavg =0.
ggavg =0.
eqwtg =0.
ymc =0.
tymc =0.
line =" "
oline=" "
C Get the name of the input file, in GHCN format. The file
C must be sorted by year (since we sum all data by month
C within a year.) The name of the output file will be that of
C the input_file.yrs.GAT where GAT stands for Global Average
C Temperature.
C2345*7891 2 3 4 5 6 712sssssss8
call getarg(1,line)
oline=trim(line)//".yrs.GAT"
open(1,file=line,form='formatted')
open(10,file=oline,form='formatted') ! output
C Read in a line of data (Country Code, ID, year, temperatures)
C Set the max year so far to this first year, set the "LASTID"
C to zero so it will fail the equality test later.
read(1,'(i3,i8,1x,i4,12i5)',end=200) icc,id,iyr,itmp
iyrmax = iyr
LASTID = 0
rewind 1
20 CONTINUE
read(1,'(i3,i8,1x,i4,12i5)',end=200) icc,id,iyr,itmp
if(iyr .gt. iyrmax) then
C if you have a new year value, you come into this loop,
C calculate the Monthly Global Average Temperatures, the
C Yearly GAT for iyrmax.
C Print it all out, and move on.
do m=1,12
if (incount(m) .ne. 0) then
C We keep a running total of tenths of degree C in itmptot,
C by month. Then we divide this by the integer count of
C valid records that went into each month. This truncates
C the result (I think this is valid, since we want to know
C conservatively how much GIStemp warmed the data
C not how much my math in this diagnostic warms the data ;-)
C So we have a "loss" of any precision beyond the "INTEGER"
C values being divided, but since they are in 1/10C, we are
C tossing 1/100C of False Precision, and nothing more.
C THEN we divide by 10. (REAL) and yield a temperature
C average for that month for that year (REAL).
C I could do a 'nint' instead: nint(itmptot(m)/incount(m))
C and get a rounded result rather than truncated, but I
C doubt if it's really worth if for a "hand tool" that I'd
C like to be a conservative one. If I truncate, then any
C "warming" of the data is from GIStemp, not this tool.
C (Or GHCN, now that I'm using to analyse the input data
C as well as the code itself.)
C2345*7891 2 3 4 5 6 712sssssss8
C Diagnostic write to check missing data flag handling.
C write(*,*) "tmpavg: ", tmpavg
tmpavg(m) = (itmptot(m)/incount(m))/10.
C Diagnostic write to check missing data flag handling.
C write(*,*) "TMPavg: ", tmpavg
gavg = gavg+tmpavg(m)
ymc = ymc+1.
C We put a running total of yearly averages together,
C along with a count for tmpavg, it is the total of
C monthly temperature averages divided by the count of
C months with data in them (converted to C from 1/10 C).
C For eqwt it is a running total of those averages that are
C used at the end to calculate a "monthly average of
C monthly averages ".
C Basically, the first form, gavg, weights each recored
C equally, while the second form gives equal weight to
C each month, regardless of number of records in that month.
C Which one is right? You get to choose... (And THAT
C is just one of the issues with an "average of averages of
C averages" means something...
C I just put them here so you can see that they are, in fact,
C different...
eqwt(m) = eqwt(m)+tmpavg(m)
eqwtc(m) = eqwtc(m)+1
end if
end do
gavg=gavg/ymc
C2345*7891 2 3 4 5 6 712sssssss8
C Write out the Year, the averages, the grand avg, and the
C number of thermometers in the year
write(10,'(i4,12f5.1,f5.1,i4)') iyrmax,tmpavg,gavg,iyc
C Diagnostic "writes", should you wish to use them.
C write(*,*) "iyc: ", iyc
C2345*7891 2 3 4 5 6 712sssssss8
C write(*,'("GAT/year: "i4,12f7.2,f7.2,i6,f7.2)') iyrmax,
C *tmpavg,gavg,iyc,ymc
C probably paranoia, but we re-zero the monthly arrays of data.
C ande pack tmpavg with missing data flags of -99
C
do m=1,12
incount(m) =0
itmptot(m) =0
tmpavg(m) =-99.
ymc =0.
end do
gavg =0.
iyc =0
LASTID =0
iyrmax =iyr
C hang on to the present year value and ...
end if
C End of "new year" record handling.
C So we have a new record (for either a new year or for the
C same year.) If it is valid data (not a missing data flag)
C add it to the running totals and increase the valid
C data count by one.
C2345*7891 2 3 4 5 6 712sssssss8
C Increment the running total for stations in this year.
C In Year Counter
if (id .NE. LASTID) then
iyc = iyc+1
LASTID = id
end if
C For each month, skipping missing data flags, increment
C the valid data counter for that month incount, add that
C temperature data (in 1/10 C as an integer) into the
C yearly running total itmptot.
C Also do the same for the total records count nncount
C and running total of all temperatures (by month) ntmptot.
do m = 1,12
if (itmp(m) .gt. -9000) then
incount(m) = incount(m)+1
itmptot(m) = itmptot(m)+itmp(m)
nncount(m) = nncount(m)+1
ntmptot(m) = ntmptot(m)+itmp(m)
end if
end do
C and go get another record
goto 20
C UNTIL we are at the end of the file.
200 continue
C2345*7891 2 3 4 5 6 712sssssss8
C Here we use the method vetted in the earlier program
C totghcn.f where we hold the temps as integers in 1/10 C
C until the very end, then we do a convert to real
C (via divide by 10.) and cast into a real (ttmpavg)
C that is the total average temperature for that month for
C the total data.
C ggavg is the grand total GAT, but after it it stuffed with
C valid data we must divide it by the number of months with
C valid data. It is the "Average of yearly averages of
C monthly averages of daily MIN/MAX averages".
C Why? Heck, GIStemp is a "serial averager", thought it might
C be fun to see what you get.
C We also show the average of all the individual monthly
C data. That gives a different value.
C Will the real GAT please stand up? ...
C I would chose to use the average of all data in a month
C since it is less sensitive to the variation of number
C of thermometers in any given year, but you might chose a
C different GAT. Averaging the data directly gives weight
C to the years with more data. Averaging the monthly
C averages gives each month equal weight. Choose one...
C Rational? No.
C But it is the reality on the ground..
C Basically, if you do "serial averaging", the order of the
C averaging will change your results. As near as I can tell,
C GIStemp (and the whole AGW movment) pay no attention to this
C Inconvenient Fact.
do m=1,12
if (nncount(m).ne.0) then
ttmpavg(m)=(ntmptot(m)/nncount(m))/10.
ggavg=ggavg+ttmpavg(m)
tymc=tymc+1
end if
eqwt(m)=eqwt(m)/eqwtc(m)
eqwtg=eqwtg+eqwt(m)
end do
C ggavg is the grand average of monthly averages for a month.
C eqwt is the sum of all months averages. eqwtc is the
C count of all months with valid data. So this is the
C place where the total gets divided by the count to give the
C average of all averages in a month.
ggavg=ggavg/tymc
eqwtg=eqwtg/tymc
write(10,'(4x,12f5.1,f5.1)') ttmpavg,ggavg
write(10,'(4x,12f5.1,f5.1)') eqwt,eqwtg
stop
end
[chiefio@tubularbells analysis]$
The wrapper script for it is:
[chiefio@tubularbells analysis]$ cat dotemps
# First off, sort v2.mean into a version for reporting by year.
DIR=${2-./Temps}
echo " "
echo -n "Do the extract / process for v2.mean_comb for ${1-501} (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
PAT=^${1-501}
echo $PAT
grep $PAT /gnuit/GIStemp/STEP0/to_next_step/v2.mean_comb > $DIR/v2.meanC.${1-501}
ls -l $DIR/v2.meanC.${1-501}
echo Now Sort
echo
sort -n -k1.13,1.16 -k1.1,1.12 $DIR/v2.meanC.${1-501} > $DIR/Temps.${1-501}
echo
echo After the Sort
echo
fi
ls -l $DIR/Temps.${1-501}
echo " "
echo "Doing GAT Yearlies w/ Missing Flag: lmyears"
echo " "
echo " "
echo -n "Do the Reporting process for $DIR/Temps.${1-501} (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
# bin/yearsghcn v2.meanC.sorted
# bin/locyearsghcn Temps.${1-501}
bin/lmyears $DIR/Temps.${1-501}
echo " " >> $DIR/Temps.${1-501}.yrs.GAT
echo For Country Code ${1-501} >> $DIR/Temps.${1-501}.yrs.GAT
echo " "
echo "Produced:"
echo " "
ls -l $DIR/Temps.${1-501}.yrs.GAT
fi
echo " "
echo -n "Look at $DIR/Temps.${1-501}.yrs.GAT (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
cat $DIR/Temps.${1-501}.yrs.GAT
fi
echo " "
ls -l $DIR/v2.meanC.${1-501} $DIR/Temps.${1-501}
echo " "
echo -n "Clean up / Delete intermediate files (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
rm $DIR/v2.meanC.${1-501} $DIR/Temps.${1-501}
fi
[chiefio@tubularbells analysis]$
Thermometer Percentages By Latitude
Once you have this combined file, you have temperature data with descriptions attached. At that time you can do “by latitude” and “by altitude” studies on the stations, countries, etc. As this “by latitude” program demonstrates:
[chiefio@tubularbells analysis]$ cat src/latcust.f
C2345*fff1 2 3 4 5 6 712sssssss8
C
C this program sorts records into latitude bands
C Input must already be sorted by year and filtered to selected latitude
C A spcial v2.mean+v2.inv concatinated file is the source
C
C There is an input file named "BANDS" that holds 9 latitude integers. S to N
C
C Copyright (c) 2009
C
C This program is free software; you can redistribute it and/or
C modify it under the terms of the GNU General Public License as
C published by the Free Software Foundation; either version 2,
C or (at your option) any later version.
C
C This program is distributed in the hope that it will be useful,
C but WITHOUT ANY WARRANTY; without even the implied warranty of
C MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
C GNU General Public License for more details.
C
C You will notice this "ruler" periodically in the code.
C FORTRAN is position sensitive, so I use this to help me keep
C track of where the first 5 "lable numbers" can go, the 6th
C position "continuation card" character, the code positions up
C to number 72 on your "punched card" and the "card serial number"
C positions that let you sort your punched cards back into a proper
C running program if you dropped the deck. (And believe it or
C not, I used that "feature" more than once in "The Days Before
C Time And The Internet Began"...
C2345*7891 2 3 4 5 6 712sssssss8
integer itmp(12), icc, id, iyr, iyrmax, m, iyc, kyr, ky, band(9)
real latcnt(11), lattot(11)
real latitude, kount
character*128 line, oline, pline
data latcnt /0,0,0,0,0,0,0,0,0,0,0/
data lattot /0,0,0,0,0,0,0,0,0,0,0/
icc=0
id=0
iyr=0
iyc=0
iyrmax=0
kyr=0
latitude=0
kount=1.
C Believe it or not, the program may make bogus values,
C unless you do this initialization.
do m=1,11
latcnt(m)=0
lattot(m)=0
end do
C Get the name of the input file, in modified GHCN format. The file must be
C sorted by year (since we sum all data by month within a year.)
C The name of the output file will be that of the inputfile.yrs.LAT
C The input file ought to be a GHCN format file with V2.inv data
C concatenated per line.
call getarg(1,line)
oline=trim(line)//".per.LAT"
pline=trim(line)//".Dec.LAT"
open(1,file=line,form='formatted')
open(10,file=oline,form='formatted') ! output
open(12,file=pline,form='formatted') ! output
open(2,file="BANDS",form='formatted') ! output
read(2,'(9i4)',end=300) band
C2345*fff1 2 3 4 5 6 712sssssss8
C Read in a line of data (Country Code, ID, year, temperatures, latitude)
C For each year, total "thermometer counts" increment decade counts
C increment thermometer counts by latitude.
write(10,'(" ")')
write(10,'(" Year SP",i4,2x,i4,2x,i4,2x,i4,2x,i4,2x,i4,2x,
*i4,2x,i4,2x,i4,"NP")') band
write(12,'(" Year SP",i4,2x,i4,2x,i4,2x,i4,2x,i4,2x,i4,2x,
*i4,2x,i4,2x,i4,"NP")') band
write(10,'(" Year SP-45 50 55 60 65 70
* 75 80 85 -NP ")')
write(12,'(" Year SP-45 50 55 60 65 70
* 75 80 85 -NP ")')
read(1,'(12x,i4,12i5,33x,f6.2)',end=200) iyr,itmp,latitude
iyrmax=iyr
rewind 1
20 read(1,'(12x,i4,12i5,33x,f6.2)',end=200) iyr,itmp,latitude
if(iyr.gt.iyrmax) then
C if you have a new year value, you come into this loop, print the
C the total thermomether count per latitude for that year
C calculate the decade total thermometer count, and every decade
C print it all out, and move on.
C2345*fff1 2 3 4 5 6 712sssssss8
C increment the latitude totals for the decade
do m=1,11
lattot(m)=lattot(m)+latcnt(m)
end do
do m=1,11
latcnt(m)=(latcnt(m)/latcnt(11))*100.
end do
C write(10,'("LAT year: "i4,9i6,1x,i5)') iyrmax,latcnt, iyc
write(10,'("LAT pct: "i4,11f6.1,1x)') iyrmax,latcnt
if (mod(iyr,10).eq.0) then
C ok, at this point we want to print out the decade average of thermometer
C counts by latitude band. In 5 degree increments.
C we would do that by printing out
kyr=iyrmax
do m=1,11
lattot(m)=((lattot(m))/lattot(11))*100.
end do
kyc=nint(kyc/kount)
write(10,'(" ")')
C write(10,'(" ",f6.1)') kount
C write(10,'("DecadeLat: "i4,9i6,1x,i5)') kyr,lattot, kyc
write(10,'("DecLatPct: "i4,11f6.1)') kyr,lattot
write(10,'(" ")')
write(12,'("DecLatPct: "i4,11f6.1)') kyr,lattot
C Then we set the decade counter to zero and reset the decade array.
do m=1,11
lattot(m)=0
end do
kount=0.
kyc=0
end if
C we re-zero the array of latitude counts for the year.
C
do m=1,11
latcnt(m)=0
end do
kount=kount+1.
iyrmax=iyr
iyc=0
end if
C2345*fff1 2 3 4 5 6 712sssssss8
C So we have a new record for a new year or for the same year.
C we count the thermometer regardless of the data flag (not many all zero)
C and we add a count to that thermometers latitude for that year.
iyc=iyc+1
kyc=kyc+1
if (latitude .lt. band(1) ) then
latcnt(1)=latcnt(1)+1
else if(latitude .lt. band(2) .and. latitude .ge. band(1) ) then
latcnt(2)=latcnt(2)+1
else if(latitude .lt. band(3) .and. latitude .ge. band(2) ) then
latcnt(3)=latcnt(3)+1
else if(latitude .lt. band(4) .and. latitude .ge. band(3) ) then
latcnt(4)=latcnt(4)+1
else if(latitude .lt. band(5) .and. latitude .ge. band(4) ) then
latcnt(5)=latcnt(5)+1
else if(latitude .lt. band(6) .and. latitude .ge. band(5) ) then
latcnt(6)=latcnt(6)+1
else if(latitude .lt. band(7) .and. latitude .ge. band(6) ) then
latcnt(7)=latcnt(7)+1
else if(latitude .lt. band(8) .and. latitude .ge. band(7) ) then
latcnt(8)=latcnt(8)+1
else if(latitude .lt. band(9) .and. latitude .ge. band(8) ) then
latcnt(9)=latcnt(9)+1
else if(latitude .ge. band(9) ) then
latcnt(10)=latcnt(10)+1
else
write(*,*) "You can't get here, compiler error Or dirty Data! "
end if
latcnt(11)=latcnt(11)+1
C and go get another record
goto 20
C UNTIL we are at the end of the file where we print the last average
200 continue
do m=1,11
lattot(m)=lattot(m)+latcnt(m)
end do
do m=1,11
latcnt(m)=(latcnt(m)/latcnt(11))*100.
end do
C write(10,'("LAT year: "i4,9i6,1x,i5)') iyrmax,latcnt, iyc
write(10,'("LAT pct: "i4,11f6.1)') iyrmax,latcnt
do m=1,11
lattot(m)=((lattot(m))/lattot(11))*100.
end do
kyc=nint(kyc/kount)
kyr=iyrmax
C write(10,'(" ",f6.1)') kount
C write(10,'("DecadeLat: "i4,9i6,1x,1i4)') kyr,lattot, kyc
write(10,'(" ")')
write(10,'("DecLatPct:"i4,11f6.1)') kyr,lattot
write(12,'("DecLatPct: "i4,11f6.1)') kyr,lattot
C2345*fff1 2 3 4 5 6 712sssssss8
300 continue
stop
end
[chiefio@tubularbells analysis]$
And the wrapper script that runs it and tends the environment:
[chiefio@tubularbells analysis]$ cat dolats
DIR=${3-./Lats}
echo " "
echo "Remember to update BANDS with 9 LAT bands prior to use"
echo " "
echo "Need to make a joined GHCN with v2.inv data for the ${1-403} records "
echo " "
ls -l $DIR/v2.${1-403}.withlat
echo " "
echo -n "Make the Extract of v2.inv.id.withlat (Y/N)? "
read ANS
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
ls -l ${2-./vetted/v2.inv.id.withlat}
echo " "
grep "^${1-403}" ${2-./vetted/v2.inv.id.withlat} > $DIR/v2.${1-403}.withlat
echo " "
ls -l $DIR/v2.${1-403}.withlat
fi
echo " "
echo "Then sort the Special GHCN with v2.inv by year"
echo "into a version for reporting."
echo " "
echo from $DIR/v2.${1-403}.withlat into $DIR/Therm.by.lat${1-403}
echo " "
ls -l $DIR/Therm.by.lat${1-403}
echo " "
echo -n "Re-sort the selected records back into year order (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
sort -n -k1.13,1.16 $DIR/v2.${1-403}.withlat > $DIR/Therm.by.lat${1-403}
ls -l $DIR/Therm.by.lat${1-403}
fi
echo " "
echo -n "Do the Count of therm/yrs by latatitude (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
echo bin/latcust $DIR/Therm.by.lat${1-403}
bin/latcust $DIR/Therm.by.lat${1-403}
echo " " >> $DIR/Therm.by.lat${1-403}.Dec.LAT
echo " " >> $DIR/Therm.by.lat${1-403}.per.LAT
echo For COUNTRY CODE: ${1-403} >> $DIR/Therm.by.lat${1-403}.Dec.LAT
echo For COUNTRY CODE: ${1-403} >> $DIR/Therm.by.lat${1-403}.per.LAT
fi
ls -l $DIR/Therm.by.lat${1-403}.Dec.LAT $DIR/Therm.by.lat${1-403}.per.LAT
echo " "
echo -n "Look at $DIR/Therm.by.lat${1-403}.Dec.LAT (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
cat $DIR/Therm.by.lat${1-403}.Dec.LAT
fi
echo " "
echo -n "Look at $DIR/Therm.by.lat${1-403}.per.LAT (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
cat $DIR/Therm.by.lat${1-403}.per.LAT
fi
echo " "
ls -l $DIR/Therm.by.lat${1-403}.Dec.LAT $DIR/Therm.by.lat${1-403}.per.LAT
echo " "
echo -n "Clean Up / Remove REPORT files (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
echo rm $DIR/Therm.by.lat${1-403}.Dec.LAT $DIR/Therm.by.lat${1-403}.per.LAT
rm $DIR/Therm.by.lat${1-403}.Dec.LAT $DIR/Therm.by.lat${1-403}.per.LAT
fi
echo " "
ls -l $DIR/v2.${1-403}.withlat $DIR/Therm.by.lat${1-403}
echo " "
echo -n "Clean Up / Remove intermediate WORK files (Y/N)? "
read ANS
echo " "
if [ "$ANS" = "Y" -o "$ANS" = "y" ]
then
echo rm $DIR/v2.${1-403}.withlat $DIR/Therm.by.lat${1-403}
rm $DIR/v2.${1-403}.withlat $DIR/Therm.by.lat${1-403}
fi
exit
echo " "
[chiefio@tubularbells analysis]$
Some of the postings about particular places and groups of countries depend on variations on these themes that select out individual countries or do a specific list. They are fairly easily created from this code base, so I’m not including them here at this time. If there is interest, I can post them too. I may do it anyway when I get time. For example, the ‘by altitude’ program is mostly a change of any LAT or lat to ALT or alt and one change of the format field to pick up data from the altitude field instead of the temperature field. Oh, and I made it an integer instead of a float. All in all, just few minutes work. I intend to merge the ALT and LAT versions into one program with a flag rather than keep two almost identical programs to maintain… so I have not posted the “by altitude” variant here, yet.
Some of the listings extend a bit past the right margin. Viewing the page source ought to let you see those bits if you need them. As time permits, I’ll come back and “pretty print” the listings so you can see the bits off the edge…
