pro csiro_annual_climat

; Purpose
; -------
; Reads monthly-mean output from the CSIRO AGCM, and generates a netCDF file
; containing a timeseries of annual means.
;
; History
; -------
; 2002 Jul 15	Steven Phipps	Original version
; 2004 Apr 14	Steven Phipps	Modified for 50-year netCDF files
; 2004 Dec 28	Steven Phipps	Tidied up, and modified for five-digit year
;				numbers
; 2005 Jan 18	Steven Phipps	Modified so as not to read the year numbers
;				from the model output, which is unnecessary
; 2009 Apr 24	Steven Phipps	Modified for the conversion of AGCM output from
;				16-bit INTEGER to 32-bit REAL

; Define parameters
nx = 64
ny = 56
nmonth = 12
nyear = 50
days = [31.0, 28.0, 31.0, 30.0, 31.0, 30.0, 31.0, 31.0, 30.0, 31.0, 30.0, 31.0]
missing_value = -1.0e34

; Get various information from user
var = ""
run_name = ""
year1 = 0
year2 = 0
read, var, prompt="Enter the variable name [e.g. tsu]  "
read, run_name, prompt="Enter the run name      [e.g. d00]  "
read, year1, prompt="Enter the initial year    [e.g. 1]  "
read, year2, prompt="Enter the final year    [e.g. 800]  "

; Convert year numbers to five-character strings
year1s = strtrim(string(year1), 2)
year2s = strtrim(string(year2), 2)
if (year1 lt 10000) then year1s = '0' + year1s
if (year1 lt 1000) then year1s = '0' + year1s
if (year1 lt 100) then year1s = '0' + year1s
if (year1 lt 10) then year1s = '0' + year1s
if (year2 lt 10000) then year2s = '0' + year2s
if (year2 lt 1000) then year2s = '0' + year2s
if (year2 lt 100) then year2s = '0' + year2s
if (year2 lt 10) then year2s = '0' + year2s

; Calculate total number of years
num_years = year2 - year1 + 1

; Calculate number of files from which to read
nfiles = fix(float(year2 - year1) / float(nyear)) + 1

; Declare arrays to hold data
clim = fltarr(nx, ny, num_years)

; Initialise output array
clim(*, *, *) = 0.0

; Loop over files
for file = 0, nfiles-1 do begin

  ; Calculate years covered by this file
  y1 = year1 + file*nyear
  y2 = y1 + nyear - 1

  ; Convert year numbers to five-character strings
  y1s = strtrim(string(y1), 2)
  y2s = strtrim(string(y2), 2)
  if (y1 lt 10000) then y1s = '0' + y1s
  if (y1 lt 1000) then y1s = '0' + y1s
  if (y1 lt 100) then y1s = '0' + y1s
  if (y1 lt 10) then y1s = '0' + y1s
  if (y2 lt 10000) then y2s = '0' + y2s
  if (y2 lt 1000) then y2s = '0' + y2s
  if (y2 lt 100) then y2s = '0' + y2s
  if (y2 lt 10) then y2s = '0' + y2s

  ; Open netCDF file
  filename = "s" + var + "_" + run_name + "_" + y1s + "-" + y2s + ".nc"
  print, "Reading data from file ", filename, " ..."
  ncid = ncdf_open(filename)

  ; Read data
  datid = ncdf_varid(ncid, var)
  ncdf_varget, ncid, datid, data

  ; If this is the first file, read the axis data and attributes
  if (file eq 0) then begin
    lonid = ncdf_varid(ncid, "longitude")
    latid = ncdf_varid(ncid, "latitude")
    ncdf_varget, ncid, lonid, lon
    ncdf_varget, ncid, latid, lat
    ncdf_attget, ncid, datid, "long_name", long_name_in
    long_name = string(long_name_in)
    vardata = ncdf_varinq(ncid, datid)
    num_atts = vardata.natts
    has_units = 0
    for i = 0, num_atts-1 do begin
      attname = ncdf_attname(ncid, datid, i)
      if (attname eq "units") then has_units = 1
    endfor
    if (has_units eq 1) then begin
      ncdf_attget, ncid, datid, "units", units_in
      units = string(units_in)
      print, "Variable ", var, " has units ", units, "."
    endif else begin
      print, "Variable ", var, " has no units."
    endelse
  endif

  ; Close netCDF file
  ncdf_close, ncid

  ; Fix missing values
  index = where(data gt 9.0e36)
  if (index ne -1) then data(index) = missing_value

  ; Derive annual means
  for month = 0, nmonth-1 do begin
    clim(*, *, y1-year1:y2-year1) = clim(*, *, y1-year1:y2-year1) + $
                                    days(month) * data(*, *, month, *) / 365.0
  endfor

endfor

; Fix missing values
index = where(clim lt -1.0e30)
if (index ne -1) then clim(index) = missing_value

; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
filename = "s" + var + "_" + run_name + "_ann_" + year1s + "-" + year2s + ".nc"
ncid = ncdf_create(filename)

; (2) Define dimensions
londid = ncdf_dimdef(ncid, "longitude", nx)
latdid = ncdf_dimdef(ncid, "latitude", ny)
yrdid = ncdf_dimdef(ncid, "year", num_years)

; (3) Define variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
yrvid = ncdf_vardef(ncid, "year", yrdid, /long)
climid = ncdf_vardef(ncid, var, [londid, latdid, yrdid], /float)

; (4) Put attributes
ncdf_attput, ncid, climid, "long_name", "Annual-mean " + long_name
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, yrvid, "units", "years"
if (has_units eq 1) then begin
  ncdf_attput, ncid, climid, "units", units
endif
ncdf_attput, ncid, climid, "missing_value", missing_value
ncdf_attput, ncid, lonvid, "modulo", " "
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, yrvid, "point_spacing", "even"
ncdf_attput, ncid, /global, "title", "CSIRO model run " + run_name + $
  ", years " + year1s + " to " + year2s

; (5) Exit define mode
ncdf_control, ncid, /endef

; (6) Write data
ncdf_varput, ncid, lonvid, lon
ncdf_varput, ncid, latvid, lat
ncdf_varput, ncid, yrvid, indgen(num_years) + year1
ncdf_varput, ncid, climid, clim

; (7) Close file
ncdf_close, ncid

end