pro csiro_monthly_climat

; Purpose
; -------
; Reads monthly-mean data from multiple CSIRO AGCM output files, and writes the
; data to a single netCDF file.
;
; History
; -------
; 2005 Mar 8	Steven Phipps	Original version
; 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
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 array to hold data
clim = fltarr(nx, ny, nmonth, num_years)

; 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

  ; Transfer data to output array
  clim(*, *, *, y1-year1:y2-year1) = data

endfor

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

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

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

; (4) Put attributes
ncdf_attput, ncid, climid, "long_name", long_name
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, monvid, "units", "months"
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, monvid, "modulo", " "
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, monvid, "point_spacing", "even"
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, monvid, indgen(nmonth) + 1
ncdf_varput, ncid, yrvid, indgen(num_years) + year1
ncdf_varput, ncid, climid, clim

; (7) Close file
ncdf_close, ncid

end