GHCN v1 vs v3 some code

These are the computer programs I use to create the v1 vs v3 comparisons.

If you are not a “computer person”, it will likely have no interest for you at all.

If you are a “computer person”, some of you will be Old Farts who learned FORTRAN decades ago and this will find a soft spot for you. Others will be Young Turks who think that only things written in Object Oriented Languages with all sorts of levels of libraries and objects calling objects can possibly be of interest. Everything procedural and especially those with things like fixed format statements were discarded decades ago. You will be wrong, but that will not trouble you in the slightest… ( There are still huge blocks of code written in COBOL, FORTRAN, Pascal, Ada, PL/1, and a few dozen other languages of that sort. Lots of data is still floating around in fixed format files, and while I can write C to handle such files, it’s a pain.)

Why FORTRAN? Because GIStemp was written in it and I was porting GIStemp when I first started working with the GHCN files. Now I’ve got a closet full of programs to do things with the data and it’s just easier to use the existing code base an hack it to a slightly new form. Yes, I write C. (And a few different database products and at one time or another I’ve either written at least one functioning program or worked with someone putting database calls into a program, in each of the languages listed above. I was a DBA Consultant helping the programming staff get into the databases using their various languages of choice. Oh, and ALGOL too. I really liked ALGOL…) At any rate, I “got over it” about “My favorite language is better than your favorite language” about 30 years ago. They all work. They all do some things better than others. While I’ve not personally written in C++ or similar object oriented languages, I did manage a development team using it (that produced 3 software patents along) and I’ve learned enough to know how it works. Still not impressed all that much. “Yeah, it works too. So?” is about my opinion. You end up with giant blocks of executable binary code, measured in Megabytes for even simple things. For this code, the largest binary is 678 K Bytes. Oh, and it’s not going to suddenly break when some dynamic library changes… once it’s compiled, it’s done. I find that comforting… Yet if a project were already written in the newest Language Du Jour, I’d just us it, whatever it was.

At any rate, with the vain hope that will hold off some of the “Language Bigotry Wars, part 237″… and knowing it can’t… I’m going to put some of the programs her for folks to see.

Yes, much of this would be easier with the data in a database. Doing all the work to pick one, install it, etc. is not high on my “to do” list right now. I had to use FORTAN to make GIStemp go, and this is now working. So why take on a load of infrastructure building that doesn’t have payback? It’s on the Someday Maybe list, but don’t hold your breath…

A couple of things, first, though.

Stealing Characters

WordPress tries to turn things into HTML. That means things like an open bracket < tends to be "stolen" and sometimes all the stuff after it. I may not catch them all on first try.

Rehashed Rewritten Reused code

This is rehashed code. It started live in 2009 processing Temperatures, not Anomalies. It’s been through the conversion and rehash wringer a couple of times, and had a dozen variations made for various special cases. (Yes, Object Language Bigots, I know you can reuse modules from your object library only writing changed layers to wrap around them… then again, who has memorized your whole catalog of objects and all their parameters? IF I really cared, and I don’t, I could make chunks of this into subroutines and call them, too. But since most of this is way less than 100 lines of active code, why bother?) So when reading the code, I may not have changed all the references to calculating a GAT Global Average Temperature to an Anomaly. In many cases variables will have names like itmp instead of itmpanomaly or ianom. The risk of breaking things from changing variable names just wasn’t worth it. That I could just change the input file to anomalies and reuse the code was more important. Just don’t be surprised if some bits are not as clean as my usual stuff.

At least I’m reusing modules and not writing everything from scratch.

Nested Structure

Long ago I got comfortable with doing some processes with UNIX script tools and some in the target languages. This also gives nice “break points” where you can write a re-usable chunk that has a few different back ends to do different things. This means that some of the ‘flow of control’ will seem odd. ( More likely “What the Hell is he doing THAT for?! It could be done in line right here more easily…”) There are two major reasons. Some of it will be legacy “cruft” that just builds up in any code base that’s been in use for a couple of years. More of it builds in flexible points where I can just change one part and make a whole new process. So, for example, I have done various “by latitude” or “by altitude” bits of code. Most of the program chunks that the same, and I just stick in a ‘if altitude is GT FOO’ line, for example.

This does mean that some things are in Shell Scripts. That means that folks can be tempted to argue the merits of C Shell, vs Borne Shell, vs BASH, vs Korn Shell, vs,… Frankly, I’m not even sure in what shell the scripts are written. I think it’s BASH as this is a Linux box; but long ago I learned to write more “generic” code that doesn’t use obscure special features and tends to run everywhere. (Except C Shell that’s “different”… so I usually avoid it unless doing maintenance on preexisting code). Yes, there will be a dozen better ways to do some of the things I do; but few of them will be simpler or more likely to just work anywhere. Even on fresh ports of “whatever” with 1/2 the features missing… It’s how I write. If it bothers you, get over it.

Wordy Comments and Pedantic Style

Having gone back to some arcane and tricky bit of code and cursed the Evil Inconsiderate Bastard who wrote such complicated and obscure stuff and left comments like “Now we do the work”… and discovered on reading the doc.block it was me… I now leave to myself copious notes in the code. I also avoid “neat tricks” that may not be remembered 4 years and 2 languages later… If that bothers you, then copy the code and edit out all the comments. They start with “C” in FORTRAN or “#” in the shell scripts and most of the folks will likely find them for more interesting that the actual code.

The Code

I concatenate the Inventory File with the Temperature Data File. This is horribly wasteful of disk space, but makes it easier for me to write some kinds of programs later. So, for example, if I want all things in certain latitude or longitude bands, I can just pick those fields out of the current record. No need to match two files on key fields. No need to code database calls. As the disk in the computer was free about a decade ago and still isn’t full, it’s not like the cost matters. It was fast, and it is convenient.

So the first bit of code just takes the INV file and matches it to the v3.mean temperature data. I named it mkinvmean for “Make Inventory Mean combined file”. In v3 the format of the temperature data changed from v2, so I used these intermediate steps to convert some of it closer to the v2 form so reuse of prior code base would be easier. This is in Shell Script form and has a lot of things where you can pass in parameters to override defaults. So, for example, the first parameter is ‘where the inv file lives’ and the second is where the v3.mean file lives. The defaults are what I use most of the time, but I can override with a simple passed parameter to try something else. For those not familiar with Unix / Linux shell scripts, they are wonderful tools that are the glue holding all the other bits together. Threaded interpreted languages with an elegance in their own right. Once written, one script can call another just by speaking its name… You will see that in the code. In the land of *Nix the symbol “..” means “one directory up the naming tree”. An executable bit of code is often stored in ./bin or ../bin (the present directory with a subdirectory of bin or ‘go up one level’ and find a directory named ‘bin’.

echo " "
echo "Optional:  sort the v3.inv type file by Country Code / Station ID"
echo " "

DIR="../data" 
DTDIR="./data"
ls -l ${1-$DIR/v3.inv} $DTDIR/v3.inv.ccid

echo " "
echo -n "Do the sort of $DIR/v3.inv data into $DTDIR/v3.inv.ccid (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     INVFILE=$DTDIR/v3.inv.ccid
     echo INVFILE= $INVFILE

     sort -n -k1.1,1.12 ${1-$DIR/v3.inv} > $INVFILE
else
     INVFILE=${1-$DIR/v3.inv} 
     echo INVFILE= $INVFILE
fi

echo " "
echo "v3.inv data from:"
echo " "
ls -l $INVFILE
echo " "

ls -l ${2-$DIR/v3.mean} $DTDIR/v3.sort.ccid

echo " "
echo "Optional:  Sort the v3.mean type file by Country code / station ID"
echo " "
echo -n "Do the sort of v3.mean into v3.sort.ccid (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     MEANFILE=$DTDIR/v3.sort.ccid
     sort -n -k1.1,1.10 ${2-$DIR/v3.mean} > $MEANFILE
else
     MEANFILE=${2-$DIR/v3.mean} 
fi

echo " "
echo "v3.mean data from: "
echo " "
ls -l $MEANFILE
echo " "

echo "Then feed the v3.inv data, sorted by cc/station ID, and "
echo "v3.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: $DTDIR/v3.mean.inv (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
      ./bin/ccaddlatv3 $INVFILE $MEANFILE
fi

echo " "
echo "Produced the list of v3.mean.inv records "
echo "(temps, with station data, sorted by CCStationID)"
echo " "

ls -l v3.mean.inv  $DTDIR/v3.mean.inv

echo " "
     ls -l $DTDIR/v3.sort.ccid $DTDIR/v3.inv.ccid
echo " "

echo -n "Remove the work files v3.sort.ccid and v3.inv.ccid (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     echo rm $DTDIR/v3.sort.ccid $DTDIR/v3.inv.ccid
     rm $DTDIR/v3.sort.ccid $DTDIR/v3.inv.ccid
fi

echo " "
echo -n "Does v3.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.4,1.11  v2.sort.id.withlat
#     meansortidyrv3 sorts on the first 15 chars.  cc,Stn,year

        echo ./bin/meansortidyrv3 v3.mean.inv giving v3.mean.inv.ccsidyr
        ./bin/meansortidyrv3 $DTDIR/v3.mean.inv $DTDIR/v3.mean.inv.ccsidyr
	echo
        ls -l $DTDIR/v3.mean.inv.ccsidyr
fi

There are a lot of meanings coded into the file names. Why? Just for me. To mark things from one program as different from another. There isn’t a great wisdom to it, nor have I made a dictionary. But, for example, v3.mean.inv.ccsidyr is a file containing data from the v3 version of GHCN, with both mean temperature based data (as opposed to MIN or MAX based) and with Country Code, Station IS, and Year as key fields. At the end of the script above, in the last block, you can see where I commented out an old “Sort” call using a “#” and substituted a small script in .bin/meansortidyrv3 that does a sort of a “mean” based file format on Station-ID and Year.

That whole “script” is really just a one line sort command. It could be put in the main script there, just like the one that was commented out. It is:

sort -rn -k1.1,1.15 ${1-"./data/v3.mean.inv"} > ${2-"./data/v3.mean.inv.rev"}

You can see that it has different defaults for the file names and sorts on the v3 sized fields. It gets used in several other bits of scripting, so make a ‘one line script’ out of it to avoid retyping it.

What does this script do? Mostly it hand holds you through the process of matching the two files on a common key and adding the Inventory inormation to each record. Yes, you could do it “long hand” but this also lets you manage intermediate files and lets you restart at any point. I usually make one of these whenever I’m going a manual process, just so I can avoid a lot of typing as I debug what to do. Then, when it’s proven right, it’s just very convenient to let it do all the typing and I just type “mkinvmean” and say Y or N to various questions.

The actual work of matching the two files on a key and printing out the combined record is done by the FORTRAN program called by name in the middle / bottom: ccaddlatv3

It is customized for v3 (thus the v3 on the end) from the original that ran on v2 data. Remember that C starts a comment line that doesn’t actually DO anything. You will notice some “diagnostic writes” commented out in the code. Some folks pull them, I leave them in for any future debugging need during some future port or re-purposing. I’ve pretty much been 100% successful at just looking at intermediate stages of the data to do debugging, so no, I don’t use a lot of other debugging tools. Never needed to.

C     Program: 	ccaddlat.f
C     Written: 	May 23, 2012
C     Author:  	E. M. Smith
C     Function:	This program matches v3.inv and v3.mean on Station IDs
C     sorted in order by cc stationID(8) and produces a merged v3.mean
C     format file.
C     So as to match station location info with  thermometer records.
C
C     Copyright (c) 2012
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     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     ayr     - The year as char for a given set of monthly data.
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, set to the input
C               file name.GAT
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*4  mod, ayr
      character*96 buff
      character*96 buff2

      data itmp    /0,0,0,0,0,0,0,0,0,0,0,0/
      data buff /"                                                  "/
      mod=""
      ayr="    "
      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="./data/v3.mean.inv"
  
C      oline=trim(line)//".withmean"

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 v3.mean record, then a v3.inv description record
C     original addlat.f has a bug, icc is read in 2 x to same variable.

      read(2,'(a11,a4,a4,a96)',end=500) csid, ayr, mod, buff2
C      write(*,*) "Before 20: ", icc, sid, ayr, mod, buff2
   20 read(1,'(a11,a96)',end=400) cid,buff
C      write(*,*) "After 20: ", icc," ", id, " ", buff
C      write(*,*) "at 20: ", sid," ",ayr," ",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!"
C2345*7891         2         3         4         5         6         712sssssss8

   70 DO WHILE (cid .eq. csid) 
         write(10,'(a11,a4,a4,a96,1x,a96)') csid,ayr,mod,buff2,buff
C         read(2,'(i6,1x,a8,a1)',end=300) idcount, sid, mod
         read(2,'(a11,a4,a4,a96)',end=300) csid,ayr, mod, buff2
C      write(*,*) "From 70: ", icc," ", sid," ", mod," ", buff2
C         read(1,'(a11,a96)',end=200) cid,buff
      END DO
      goto 40

   80 DO WHILE (cid .gt. csid) 
C         read(2,'(i6,1x,a8,a1)',end=300) idcount, sid, mod
C         read(2,'(i3,a8,a1,a64)',end=300) icc, sid, mod, buff2
         read(2,'(a11,a4,a4,a96)',end=300) csid, ayr, 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,a96)',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 v3.inv Station Records - ought to be rare!"

  300 continue
      STOP "Out of v3.mean records.  Most likely case."

  400 STOP "Input file 1 blank on first record!"
  500 STOP "Input file 2 blank on first record!"
      END

Pretty short and fairly straight forward. Run your “digital finger” down each file, matching lines, write out the composite, move on to the next one.

Making Anomalies

This next set just turns temperatures into anomalies. It does this with a simple variation on First Differences ( a peer reviewed method, IIRC). It normally starts with the first data item, and sets that ‘difference’ to zero (as it is what it is) then looks at the next data item in the time series and subtracts one from the other, giving the difference. So if it was 20 degrees yesterday, and 21 today, the difference is 1 degree. So you get “0”, then “1” as the entries.

I make two changes.

1) I “start time” with the present record and go backwards in time. This is still a “First Differences”, but it uses the most recent temperature as the starting point. (Though I’ve sometimes called it “Last Differences” for a shorthand term). In this way we avoid the large swings that can be caused if the first data point is a bit whacky. Very old thermometers can sometimes be significantly “off” from the present ones. So ‘error’ tends to show up in the far past where folks don’t really are. It also means more thermometer records “start time” at about the same place. Convenient both for giving more consistent results, but also for highlighting when a large group of thermometers suddenly enter (going backward in time, exiting going forward in time) the records. Finally, but starting “now” and going back in time we can clearly see if the current ‘cycle wiggle’ is biasing the results. If most of the record is 1/2 C lower than the last couple of years, your “warming” is from where the record ends in a cycle, not some long term change of temperatures in a gradual way. We’ve seen that in several of the graphs.

2) Missing data handling. For classical First Differences, if there are missing data, you just re-zero and reset. As there are a lot of missing data in some places, that gives some very strange results. I just assume that if you have January for some years in a row, and leave out a couple, then get some more, that somewhere in between you passed through the space between the two. So if you had a series that was +1/10, +1/10, -1/10, missing, missing, +2/10; I just hold onto the last anomaly and wait while skipping missing data. When I hit the +2, I just account for all the change in that year. So you would have 1/10, 2/10, 1/10, 3/10 as the series. This, IMHO, more accurately reflects the reality of the temperatures recorded. That last year WAS 2/10 higher than the prior value, so throwing it away and using a zero is just wrong. In this way this code is much more robust to data dropouts and also more accurate.

We again start off with a ‘wrapper script” to do the managing of the process. I take the default values for file names and such and stuff them into symbolic names so that my typing is reduced, the code is more readable, and they can be changed in one place without risk of typos further down in the code at a later time. We also see the reuse of that one line sorting script from above. Then a call to the FORTRAN program that does the work, passing it a correctly sorted file.

#
# Take a composite mean and inventory file, and reverse sort it.
# Then feed it to the anomaly file creating program dTdtv1
#
#


DTDIR="./data"
MEANINV=${1-$DTDIR/v3.mean.inv} 
REVD=$DTDIR/v3.mean.invR

DIR="../data"

echo " "
echo -n "Do the sort of $MEANINV data into $REVD (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     INVFILE=$MEANINV
     echo INVFILE= $INVFILE
     echo
     echo ./bin/meansortmeaninv3 $MEANINV giving $REVD
     echo

     ./bin/meansortidyrv3 $MEANINV $REVD

     echo

     ls -l $REVD
fi

./bin/dTdtv3 $REVD

ls -l $REVD.dt

And what does dTdtv3 look like? It has a bit fewer comments, but is very short. I’ve also got some ‘commented out’ diagnostic writes that print descriptive things when active, but act as document / comments when commented out. A useful “trick”.

C2345*fff1         2         3         4         5         6         712sssssss8
C
C    this program converts temperatures into delta Temp / delta time
C    Input must already be sorted by Station ID and by year 
C    A special v3.mean+v3.inv concatenated file is the source (so that
C    station meta-data are attached to temperatures
C

      integer temps(12), otemps(12), dtemps(12)
      integer yr,  m, bad

      character*11 stnID, ostnID 
      character*128 meta, oline, line

      data temps   /0,0,0,0,0,0,0,0,0,0,0,0/
      data otemps  /0,0,0,0,0,0,0,0,0,0,0,0/
      data dtemps  /0,0,0,0,0,0,0,0,0,0,0,0/

      yr     = 0 
      m      = 0
      bad    = -9999

      stnID  = " "
      ostnID = " "
      meta   = " "
      oline  = " "
      line   = " " 


C     Get the name of the input file, in modified GHCN format (w/ metadata).
C     The file must be sorted by year (since we do a delta year to year.)
C     The name of the output file will be that of the inputfile.dt
C     The input file ought to be a GHCN format file with V2.inv data
C     concatenated at the end of each temperature line, as will the 
C     output file have meta data after each dT/dt line.

      call getarg(1,line)
      oline=trim(line)//".dt"
      open(1,file=line,form='formatted')
      open(10,file=oline,form='formatted')              ! output


C2345*fff1         2         3         4         5         6         712sssssss8

C     Read in a line of data 
C     (Country Code/StnID, year, temperatures, meta)

   20 read(1,'(a11,i4,4x,12(i5,3x),1x,a96)',end=999) stnID,yr,temps,meta
C     write (*,*) "Input Data"
C     write (*,*) stnID, yr, temps, meta

C      if you have a new record for the same station, find dT/dt.

C2345*fff1         2         3         4         5         6         712sssssss8

      if(stnID .eq. ostnID) then
        do m=1,12
           if ( (otemps(m).ne.bad) .and. (temps(m).ne.bad) ) then
C              write (*,*) "non-bad temp and otemp"
               dtemps(m) = temps(m) - otemps(m)
               otemps(m) = temps(m)
           else
C              write (*,*) "Somebody Bad "
               dtemps(m)=0
               if ( temps(m) .ne. bad ) then
C                   write (*,*) "temps non-bad +++++"
                    otemps(m) = temps(m)
               end if
           end if

         end do

      else
           ostnID=stnID
           do m=1,12
              dtemps(m)=0
              otemps(m)=temps(m)
           end do
      end if 
C     write (*,*) "ORDINARY Output"
C     write (*,*) stnID, yr, dtemps, meta
      write (10,'(a11,i4,12i5,a96)') stnID, yr, dtemps, meta

      go to 20

  999 continue
C     write (*,*) "Final RECORD"
C     write (*,*) stnID, yr, dtemps, meta

      close (1)
      close (10)

C2345*fff1         2         3         4         5         6         712sssssss8
 9999 continue

      stop
      end

Again, too, notice the comment – “rulers” that give column positions. On the old punched cards, spaces in the front were reserved for line numbers. At position 72, things just ‘cut off’ (and if you don’t use a proper ‘continuation card’ can cause odd bugs when the variable “alltojohn” becomes “alltojo” as the last to characters get chopped…) So I have those ‘rulers’ stuck in every so often to make it easy to spot lines that need a continuation…

This program is compiled with it’s own little 2 line script. It is named “mkdTf”

g95 src/dTdtv3.f 
mv a.out bin/dTdtv3 

In this way I don’t have to worry about which compiler I used (f77? g95?). I can also just comment out the last line when I want to experiment with new versions without overwriting the old one that works.

There’s a similar program for aligning the v3 data with the v1 data (that ends in 1990). It differs by one line. Just after that line 20 that reads in the data, it looks at the year and if it is newer than 1990, just goes to read the next record. Here is the output of the ‘diff’ program that compares the two:

54a55
>       if(yr .gt. 1990) goto 20

Once the file is built, it’s easy to make many different kinds of reports from it. Here’s one of them. The dPdt0 program that looks at data which ends in 1990 and makes a report. First up, the script, that makes sure the right files are being used. Again, it’s an interactive hand holder that cuts down on typing.

It also does nice little coordination things like letting me choose to keep the intermediate files around (so future reports are faster) or remove them; and lets me see the report (in ‘vi’ the Visual Editor) and asks me if I want to put them on my ftp server so I can make graphs out of them on the laptop or not.

#       This version expects a file with the 'mod' flags combined.
#       such as ought to be in ./data/v3.mean.inv.dt

#	First off, extract the target country or countries data from a v3.mean type 
#       file, then sort it into a version for reporting by year.

DIR=${3-./DTemps}
DAT=${4-./data}
INFILE=${2-$DAT/v3.mean.invR0.dt}
STATN=ALL

TEMPS=Temps.rM0$STATN
REPORT=$DIR/$TEMPS.yrs.dP
#STATN=${1-501}

OUTFILE=$DIR/v3.meanR0.$STATN

echo " " 
echo "Do the source and  Work Files already exist? : "
echo " "
ls -l $INFILE  $OUTFILE $DIR/$TEMPS

echo " "
echo "Do Older Reports in this series exist? : "
echo " "
ls -l $REPORT
echo " "
echo -n "Extract / process $STATN from $INFILE (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     PAT=^[0-8]
     echo $PAT
     grep $PAT $INFILE > $OUTFILE

     ls -l $OUTFILE
fi

echo " "
echo -n "Do the SORT for $OUTFILE (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     echo Now Sort
     echo

     sort -rn -k1.12,1.15 -k1.1,1.11  $OUTFILE > $DIR/$TEMPS

     echo 
     echo After the Sort
     echo 
fi

ls -l $DIR/$TEMPS

echo " "
echo "Doing dP/dt Yearlies: dpyears"
echo " "

echo " "
echo -n "Do the Reporting from $DIR/$TEMPS (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
#     ./bin/lrmdyearsv3 $DIR/$TEMPS
#     ./bin/dMTdt $DIR/$TEMPS
#     ./a.out $DIR/$TEMPS

     ./bin/dpyears $DIR/$TEMPS

     echo " " >> $REPORT
     echo For Country Code $STATN >> $REPORT
     echo " " >> $REPORT
     echo From input file $INFILE  >> $REPORT
     echo " "
     echo "Produced:"
     echo " "
     ls -l $REPORT
fi

echo " "
echo -n "Look at $REPORT (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
#     cat $REPORT
     vi $REPORT
fi

echo " "
     ls -l $OUTFILE $DIR/$TEMPS
echo " "
echo -n "Clean up / Delete intermediate files (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     rm $OUTFILE $DIR/$TEMPS
fi

#CoStns ${1-501}

echo "Do form country name"

echo " "
echo -n "Copy $REPORT to /var/ftp/pub/ToPc/v3.$TEMPS.yrs.dP (Y/N)? "
read ANS
echo " "

if [ "$ANS" = "Y" -o "$ANS" = "y" ] 
then
     cp $REPORT /var/ftp/pub/ToPc/v3.$TEMPS.yrs.dP
     cp fort.11 /var/ftp/pub/ToPc/v3.$TEMPS.yrs.dP.csv
     ls -l  $REPORT /var/ftp/pub/ToPc/v3.$TEMPS.yrs.dP
     ls -l  /var/ftp/pub/ToPc/v3.$TEMPS.yrs.dP.csv
fi

Notice that I produce a cut down version of the summary columns in “Comma Separated Value” or csv form. Makes it a little quicker to do the graphs. It goes to the default file name “fort.11” and I rename it something more useful.

You can also see where this same “wrapper” with minor changes has been used for several other programs that are commented out with a “#” sign. This one calls ‘dpyears’ to do the actual report. In some ways this is the ugliest of the batch. It started life as a program that averages temperatures, so it still is commented with “temperatures’ in stead of anomalies. That needs polishing…

C2345*7891         2         3         4         5         6         712sssssss8
C     Program:  dpyears.f
C     Written:  May 24, 2012
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) 2012
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 FORTRAN
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, one 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     mod     - "modification history flag".  Place holder, not used.
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

C    runt     - An added running total annual anomaly... 

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, mod

      real tmpavg(12), ttmpavg(12), eqwt(12), eqwtc(12)
      real gavg, ggavg, ymc, tymc, eqwtg, runt

      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
C      mike  =0

      gavg  =0.
      ggavg =0.
      eqwtg =0.
      ymc   =0.
      tymc  =0.
      runt  =0.

      line =" "
      oline=" "

C     Get the name of the input file, in GHCN format.  The file
C     must be sorted by reverse 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.dT where dT stands for delta
C     Temperature.

C2345*7891         2         3         4         5         6         712sssssss8

      call getarg(1,line)
      oline=trim(line)//".yrs.dP"
      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.
C     Note that we read in an 8 digit station ID, and do not
C     count the individual "modification flags" as new stations.

      read(1,'(i3,i8,i4,12i5)',end=200) icc,id,iyr,itmp
      iyrmax = iyr
      LASTID = 0
      rewind 1
      write(10,'("Produced from input file: ",a)') line
      write(10,'(" ")')
      write(10,14)
      write(10,15)
      write(10,16)
      write(10,17)
      write(10,18)

   14 FORMAT("Thermometer Records, Average of Monthly dP/dt, Yearly runn
     *ing total")
C   14 FORMAT("Thermometer Records, Average of Monthly Data and Yearly Av
C     *erage")
   15 FORMAT("by Year Across Month, with a count of thermometer records
     *in that year")
   16 FORMAT("----------------------------------------------------------
     *-------------------------")
C2345*7891         2         3         4         5         6         712sssssss8
   17 FORMAT("YEAR     dP dT/yr  Count JAN  FEB  MAR  APR  MAY  JUN JULY 
     *  AUG SEPT  OCT  NOV  DEC")
   18 FORMAT("----------------------------------------------------------
     *-------------------------")

   20 CONTINUE

      read(1,'(i3,i8,i4,12i5)',end=200) icc,id,iyr,itmp
C      mike=mike+1
C      write(*,*) icc,id,iyr,itmp

      if(iyr .lt. iyrmax) then

C      if you have a new year value, you come into this loop,
C      calculate the Monthly Average delta  Temperatures, the
C      Yearly dT running total 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/100C, we are
C      tossing 1/100C of False Precision, and nothing more.
C      THEN we divide by 100. (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))/100.

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/100 C).
C       For eqwt it is a running total of those averages that are
C       used at the end to calculate a "by month 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
        runt=runt+gavg

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

C        write(10,'(i4,12f5.1,f5.1,1x,i4,1x,f6.3)') iyrmax,tmpavg,gavg,
C     *iyc,runt
C2345*7891         2         3         4         5         6         712sssssss8
        write(10,'(i4,1x,f6.2,f6.2,1x,i4,1x,12f5.1)') iyrmax,runt,gavg,
     *iyc ,tmpavg
        write(11,'(i4,",",f6.2,",",f6.2,",",i4)') iyrmax,runt,gavg,iyc

C     Diagnostic "writes", should you wish to use them.
C       write(*,*) "iyc: ", iyc, "Mike: ", mike
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)  =-999.
          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, subtract that
C     year dT as we are running time backwards...
C     temperature data (in 1/100 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. -9900) 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

C     Adding last year handling

        do m=1,12

          if (incount(m) .ne. 0) then
            tmpavg(m) = (itmptot(m)/incount(m))/100.
            gavg      = gavg+tmpavg(m)
            ymc       = ymc+1.

            eqwt(m)   = eqwt(m)+tmpavg(m)
            eqwtc(m)  = eqwtc(m)+1

          end if
        end do

        gavg=gavg/ymc
        runt=runt+gavg

C        write(10,'(i4,12f5.1,f5.1,1x,i4)') iyrmax,tmpavg,gavg,iyc
C        write(10,'(i4,12f5.1,f5.1,1x,i4,1x,f6.3)') iyrmax,tmpavg,gavg,
C     *iyc,runt

        write(10,'(i4,1x,f6.2,f6.2,1x,i4,1x,12f5.1)') iyrmax,runt,gavg,
     *iyc ,tmpavg

C2345*7891         2         3         4         5         6         712sssssss8

C     At the end, we make two "grand averages" showing how the order
C     in which an average of avereages of data is done changes the 
C     value.  This effect seems to be studiously ignored in 
C     discussions of GIStemp and other climate data "serial averagers".

C     Here we use the method vetted in the earlier program
C     totghcn.f where we hold the temps as integers in 1/100 C
C     until the very end, then we do a convert to real
C     (via divide by 100.) 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 is 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 but puts too 
C     much weight on individual thermometers in months with
C     few valid thermometers (like early years).  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.  They just assume that the "anomaly 
C     process" will save them from it somehow.


      do m=1,12
          if (nncount(m).ne.0) then
            ttmpavg(m)=(ntmptot(m)/nncount(m))/100.
            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 data.  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

C     "At" is the "Total delta data"
C     "Aa" is the "Average of monthly dT"

C      write(10,'("At",2x,12f5.2,f5.2)') ttmpavg,ggavg
C      write(10,'("At",2x,12f5.2,f5.2)') ntmptot,ggavg
C      write(10,'("Aa",2x,12f5.2,f5.2)') eqwt,eqwtg

      stop
      end

In Conclusion

Well, that’s the basic code. Look it over. Let me know if there are any errors or omissions you find. If folks want any more of the variations just let me know (like the one that only goes to 2009 to compare with v2).

It does still need a bit of clean-up, but it is what it is.

As I’ve been up all night making sure a rogue cat leaves the bunnies alone, I’m very much needing a nap right now. For that reason, My QA on these listings will have to waite a few hours… for when I can see straight again ;-)

Subscribe to feed

Advertisement