pro cdepthm_annual

; Purpose
; -------
; Reads in annual means of CSIRO ocean model output, and produces a single
; output file containing the data for the maximum depth of convection
;
; History
; -------
; 2004 Mar 4	Steve Phipps	Original version
; 2004 Sep 17	Steve Phipps	Modified for five-digit year numbers

; Define parameters
nx = 64
ny = 56
missing_value = -9999.9

; Get details from user
run_name = ""
year1 = 0
year2 = 0
read, run_name, prompt="Please enter the run name [e.g. d30]    "
read, year1, prompt="Please enter the initial year [e.g. 1]  "
read, year2, prompt="Please enter the final year [e.g. 500]  "

; Get number of years of data to read
nyears = year2 - year1 + 1

; Declare arrays to hold data
cdepthm = fltarr(nx, ny, 50)

; Derive name of output files
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
outfile = "cdepthm." + run_name + "." + year1s + "-" + year2s + ".nc"

; Create output file
;
; (1) Create netCDF file
outid = ncdf_create(outfile)

; (2) Define dimensions
londid = ncdf_dimdef(outid, "longitude", nx)
latdid = ncdf_dimdef(outid, "latitude", ny)
timdid = ncdf_dimdef(outid, "year", nyears)

; (3) Define variables
lonvid = ncdf_vardef(outid, "longitude", londid, /float)
latvid = ncdf_vardef(outid, "latitude", latdid, /float)
timvid = ncdf_vardef(outid, "year", timdid, /short)
cdepthmid = ncdf_vardef(outid, "cdepthm", [londid, latdid, timdid], /float)

; (4) Put attributes
title = "CSIRO model run " + run_name + ", years " + year1s + " to " + year2s
ncdf_attput, outid, cdepthmid, "long_name", "Maximum depth of convection"
ncdf_attput, outid, lonvid, "units", "degrees_east"
ncdf_attput, outid, latvid, "units", "degrees_north"
ncdf_attput, outid, timvid, "units", "years"
ncdf_attput, outid, cdepthmid, "units", "m"
ncdf_attput, outid, cdepthmid, "missing_value", missing_value
ncdf_attput, outid, lonvid, "modulo", " "
ncdf_attput, outid, lonvid, "point_spacing", "even"
ncdf_attput, outid, latvid, "point_spacing", "uneven"
ncdf_attput, outid, timvid, "point_spacing", "even"
ncdf_attput, outid, /global, "title", title

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

; Loop over years in 50-year increments
for year = year1, year2, 50 do begin

  ; Derive name of input file
  yr1 = year
  yr2 = year + 49
  yr1s = strtrim(string(yr1), 2)
  yr2s = strtrim(string(yr2), 2)
  if (yr1 lt 10000) then yr1s = '0' + yr1s
  if (yr1 lt 1000) then yr1s = '0' + yr1s
  if (yr1 lt 100) then yr1s = '0' + yr1s
  if (yr1 lt 10) then yr1s = '0' + yr1s
  if (yr2 lt 10000) then yr2s = '0' + yr2s
  if (yr2 lt 1000) then yr2s = '0' + yr2s
  if (yr2 lt 100) then yr2s = '0' + yr2s
  if (yr2 lt 10) then yr2s = '0' + yr2s
  infile = "com.ann." + run_name + "." + yr1s + "-" + yr2s + ".nc"
  print, "Reading data from file ", infile, " ..."

  ; Open input file
  inid = ncdf_open(infile)

  ; If this is the first input file, read the axes, and write them to the
  ; output file
  if (year eq year1) then begin
    yrs = indgen(nyears) + year1
    loninid = ncdf_varid(inid, "lonts")
    latinid = ncdf_varid(inid, "latts")
    ncdf_varget, inid, loninid, lon
    ncdf_varget, inid, latinid, lat
    ncdf_varput, outid, lonvid, lon
    ncdf_varput, outid, latvid, lat
    ncdf_varput, outid, timvid, yrs
  endif

  ; Read data from input file
  cdepthminid = ncdf_varid(inid, "cdepthm")
  ncdf_varget, inid, cdepthminid, cdepthm

  ; Write data to output file
  ncdf_varput, outid, cdepthmid, cdepthm, count=[nx, ny, 50], $
    offset=[0, 0, year-year1]

endfor

; Close output file
ncdf_close, outid

end