pro csiro_annual_extent

; Purpose
; -------
; Reads monthly-mean sea ice concentrations from the CSIRO AGCM, and generates
; a netCDF file containing a timeseries of annual-mean sea ice extent.
;
; Note that, for each month, the extent is regarded as being 1.0 if the ice
; concentration is greater than or equal to 15%, and 0.0 otherwise.
;
; History
; -------
; 2005 Nov 6	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
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
run_name = ""
year1 = 0
year2 = 0
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
extent = fltarr(nx, ny, nmonth, nyear)
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 = "sico_" + run_name + "_" + y1s + "-" + y2s + ".nc"
  print, "Reading data from file ", filename, " ..."
  ncid = ncdf_open(filename)

  ; Read data
  datid = ncdf_varid(ncid, "ico")
  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
  endif

  ; Close netCDF file
  ncdf_close, ncid

  ; Derive the sea ice extent for each month
  extent(*, *, *, *) = 0.0
  for l = 0, nyear-1 do begin
    for k = 0, nmonth-1 do begin
      for j = 0, ny-1 do begin
        for i = 0, nx-1 do begin
          if (data(i, j, k, l) ge 0.15) then extent(i, j, k, l) = 1.0
        endfor
      endfor
    endfor
  endfor

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

endfor

; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
filename = "sext_" + 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, "ext", [londid, latdid, yrdid], /float)

; (4) Put attributes
ncdf_attput, ncid, climid, "long_name", "Annual-mean sea ice extent"
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, yrvid, "units", "years"
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