pro sst_monthly

; Purpose
; -------
; Reads in monthly-mean CSIRO ocean model output, and produces a single output
; containing the SST data.
;
; History
; -------
; 2005 Jul 9	Steve Phipps	Original version

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

; Declare arrays
temp_in = fltarr(nx, ny, nz, 12)

; Get details from user
run_name = ""
year1 = 0
year2 = 0
read, run_name, 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]  "

; Declare output array
num_years = year2 - year1 + 1
sst = fltarr(nx, ny, 12, num_years)

; Loop over years
for year = year1, year2 do begin

  ; Derive name of input file
  yr = year
  yrs = strtrim(string(yr), 2)
  if (yr lt 10000) then yrs = '0' + yrs
  if (yr lt 1000) then yrs = '0' + yrs
  if (yr lt 100) then yrs = '0' + yrs
  if (yr lt 10) then yrs = '0' + yrs
  infile = "com." + run_name + "." + yrs + ".nc"
  print, "Reading data from file ", infile, " ..."

  ; Open input file
  ncid = ncdf_open(infile)

  ; If this is the first year, read the axis data
  if (year eq year1) then begin
    lonid = ncdf_varid(ncid, "lonts")
    latid = ncdf_varid(ncid, "latts")
    ncdf_varget, ncid, lonid, lon
    ncdf_varget, ncid, latid, lat
  endif

  ; Read temperature data from input file
  tempid = ncdf_varid(ncid, "temp")
  ncdf_varget, ncid, tempid, temp_in

  ; Transfer SST data to output array
  sst(*, *, *, year-year1) = temp_in(*, *, 0, *)

  ; Close input file
  ncdf_close, ncid

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
outfile = "sst.mth." + run_name + "." + year1s + "-" + year2s + ".nc"
outid = ncdf_create(outfile)

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

; (3) Define variables
lonvid = ncdf_vardef(outid, "longitude", londid, /float)
latvid = ncdf_vardef(outid, "latitude", latdid, /float)
monvid = ncdf_vardef(outid, "month", mondid, /short)
yrvid = ncdf_vardef(outid, "year", yrdid, /short)
sstid = ncdf_vardef(outid, "sst", [londid, latdid, mondid, yrdid], /float)

; (4) Put attributes
ncdf_attput, outid, sstid, "long_name", "Sea surface temperature"
ncdf_attput, outid, lonvid, "units", "degrees_east"
ncdf_attput, outid, latvid, "units", "degrees_north"
ncdf_attput, outid, monvid, "units", "months"
ncdf_attput, outid, yrvid, "units", "years"
ncdf_attput, outid, sstid, "units", "degC"
ncdf_attput, outid, sstid, "missing_value", missing_value
ncdf_attput, outid, lonvid, "modulo", " "
ncdf_attput, outid, monvid, "modulo", " "
ncdf_attput, outid, lonvid, "point_spacing", "even"
ncdf_attput, outid, latvid, "point_spacing", "uneven"
ncdf_attput, outid, monvid, "point_spacing", "even"
ncdf_attput, outid, yrvid, "point_spacing", "even"
ncdf_attput, outid, /global, "title", "CSIRO model run " + run_name + $
  ", years " + year1s + " to " + year2s

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

; (6) Write data
ncdf_varput, outid, lonvid, lon
ncdf_varput, outid, latvid, lat
ncdf_varput, outid, monvid, indgen(12) + 1
ncdf_varput, outid, yrvid, indgen(num_years) + year1
ncdf_varput, outid, sstid, sst

; (7) Close output file
ncdf_close, outid

end