pro csiro_annual_maximum

; Purpose
; -------
; Reads monthly-mean output from the CSIRO AGCM, and generates a netCDF file
; containing a timeseries of annual maxima.
;
; Note that it is only particularly meaningful to use this program on the
; variable THM, which contains the extreme maximum temperature for each month.
;
; History
; -------
; 2005 Jul 30	Steven Phipps	Original version, based on
;				csiro_annual_climat.pro
; 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 = ""
year1 = 0
year2 = 0
read, var, prompt="Enter the variable name [e.g. thm]  "
read, run, prompt="Enter the run name      [e.g. d60]  "
read, year1, prompt="Enter the initial year    [e.g. 1]  "
read, year2, prompt="Enter the final year   [e.g. 1000]  "

; 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 output array
clim = fltarr(nx, ny, 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 + "_" + 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 1.0e36)
  if (index ne -1) then data(index) = missing_value

  ; Derive annual means
  for k = 0, nyear-1 do begin
    for j = 0, ny-1 do begin
      for i = 0, nx-1 do begin
        if (min(data(i, j, *, k))) gt -1.0e30 then begin
          clim(i, j, file*nyear+k) = max(data(i, j, *, k))
        endif else begin
          clim(i, j, file*nyear+k) = missing_value
        endelse
      endfor
    endfor
  endfor

endfor

; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
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
filename = "s" + var + "_" + run + "_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 " + 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 + ", 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