pro ts_annual

; Purpose
; -------
; Reads in annual means of CSIRO ocean model output, and produces a single
; output file containing the data for each of the following variables:
; (1) surface heat flux
; (2) surface salinity tendency
; (3) temperature
; (4) salinity
; (5) SST
; (6) SSS
;
; History
; -------
; 2003 Dec 12	Steve Phipps	Original version
; 2004 Sep 17	Steve Phipps	Modified for five-digit year numbers; changed
;				units of depth from "metres" to "m"
; 2005 Jun 17	Steve Phipps	Modified to save SST and SSS as well

; Define parameters
nx = 64
ny = 56
nz = 21
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
stfht = fltarr(nx, ny, 50)
stfsal = fltarr(nx, ny, 50)
temp = fltarr(nx, ny, nz, 50)
sal = fltarr(nx, ny, nz, 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
outfile1 = "stfht." + run_name + "." + year1s + "-" + year2s + ".nc"
outfile2 = "stfsal." + run_name + "." + year1s + "-" + year2s + ".nc"
outfile3 = "temp." + run_name + "." + year1s + "-" + year2s + ".nc"
outfile4 = "sal." + run_name + "." + year1s + "-" + year2s + ".nc"
outfile5 = "sst." + run_name + "." + year1s + "-" + year2s + ".nc"
outfile6 = "sss." + run_name + "." + year1s + "-" + year2s + ".nc"

; Create output files
;
; (1) Create netCDF files
outid1 = ncdf_create(outfile1)
outid2 = ncdf_create(outfile2)
outid3 = ncdf_create(outfile3)
outid4 = ncdf_create(outfile4)
outid5 = ncdf_create(outfile5)
outid6 = ncdf_create(outfile6)

; (2) Define dimensions
londid1 = ncdf_dimdef(outid1, "longitude", nx)
londid2 = ncdf_dimdef(outid2, "longitude", nx)
londid3 = ncdf_dimdef(outid3, "longitude", nx)
londid4 = ncdf_dimdef(outid4, "longitude", nx)
londid5 = ncdf_dimdef(outid5, "longitude", nx)
londid6 = ncdf_dimdef(outid6, "longitude", nx)
latdid1 = ncdf_dimdef(outid1, "latitude", ny)
latdid2 = ncdf_dimdef(outid2, "latitude", ny)
latdid3 = ncdf_dimdef(outid3, "latitude", ny)
latdid4 = ncdf_dimdef(outid4, "latitude", ny)
latdid5 = ncdf_dimdef(outid5, "latitude", ny)
latdid6 = ncdf_dimdef(outid6, "latitude", ny)
depdid3 = ncdf_dimdef(outid3, "depth", nz)
depdid4 = ncdf_dimdef(outid4, "depth", nz)
timdid1 = ncdf_dimdef(outid1, "year", nyears)
timdid2 = ncdf_dimdef(outid2, "year", nyears)
timdid3 = ncdf_dimdef(outid3, "year", nyears)
timdid4 = ncdf_dimdef(outid4, "year", nyears)
timdid5 = ncdf_dimdef(outid5, "year", nyears)
timdid6 = ncdf_dimdef(outid6, "year", nyears)

; (3) Define variables
lonvid1 = ncdf_vardef(outid1, "longitude", londid1, /float)
lonvid2 = ncdf_vardef(outid2, "longitude", londid2, /float)
lonvid3 = ncdf_vardef(outid3, "longitude", londid3, /float)
lonvid4 = ncdf_vardef(outid4, "longitude", londid4, /float)
lonvid5 = ncdf_vardef(outid5, "longitude", londid5, /float)
lonvid6 = ncdf_vardef(outid6, "longitude", londid6, /float)
latvid1 = ncdf_vardef(outid1, "latitude", latdid1, /float)
latvid2 = ncdf_vardef(outid2, "latitude", latdid2, /float)
latvid3 = ncdf_vardef(outid3, "latitude", latdid3, /float)
latvid4 = ncdf_vardef(outid4, "latitude", latdid4, /float)
latvid5 = ncdf_vardef(outid5, "latitude", latdid5, /float)
latvid6 = ncdf_vardef(outid6, "latitude", latdid6, /float)
depvid3 = ncdf_vardef(outid3, "depth", depdid3, /float)
depvid4 = ncdf_vardef(outid4, "depth", depdid4, /float)
timvid1 = ncdf_vardef(outid1, "year", timdid1, /short)
timvid2 = ncdf_vardef(outid2, "year", timdid2, /short)
timvid3 = ncdf_vardef(outid3, "year", timdid3, /short)
timvid4 = ncdf_vardef(outid4, "year", timdid4, /short)
timvid5 = ncdf_vardef(outid5, "year", timdid5, /short)
timvid6 = ncdf_vardef(outid6, "year", timdid6, /short)
stfhtid = ncdf_vardef(outid1, "stfht", [londid1, latdid1, timdid1], /float)
stfsalid = ncdf_vardef(outid2, "stfsal", [londid2, latdid2, timdid2], /float)
tempid = ncdf_vardef(outid3, "temp", [londid3, latdid3, depdid3, timdid3], $
                     /float)
salid = ncdf_vardef(outid4, "sal", [londid4, latdid4, depdid4, timdid4], $
                    /float)
sstid = ncdf_vardef(outid5, "sst", [londid5, latdid5, timdid5], /float)
sssid = ncdf_vardef(outid6, "sss", [londid6, latdid6, timdid6], /float)

; (4) Put attributes
title = "CSIRO model run " + run_name + ", years " + year1s + " to " + year2s
ncdf_attput, outid1, stfhtid, "long_name", "Surface heat flux"
ncdf_attput, outid2, stfsalid, "long_name", "Surface salinity tendency"
ncdf_attput, outid3, tempid, "long_name", "Potential temperature"
ncdf_attput, outid4, salid, "long_name", "Salinity"
ncdf_attput, outid5, sstid, "long_name", "Sea surface temperature"
ncdf_attput, outid6, sssid, "long_name", "Sea surface salinity"
ncdf_attput, outid1, lonvid1, "units", "degrees_east"
ncdf_attput, outid2, lonvid2, "units", "degrees_east"
ncdf_attput, outid3, lonvid3, "units", "degrees_east"
ncdf_attput, outid4, lonvid4, "units", "degrees_east"
ncdf_attput, outid5, lonvid5, "units", "degrees_east"
ncdf_attput, outid6, lonvid6, "units", "degrees_east"
ncdf_attput, outid1, latvid1, "units", "degrees_north"
ncdf_attput, outid2, latvid2, "units", "degrees_north"
ncdf_attput, outid3, latvid3, "units", "degrees_north"
ncdf_attput, outid4, latvid4, "units", "degrees_north"
ncdf_attput, outid5, latvid5, "units", "degrees_north"
ncdf_attput, outid6, latvid6, "units", "degrees_north"
ncdf_attput, outid3, depvid3, "units", "m"
ncdf_attput, outid4, depvid4, "units", "m"
ncdf_attput, outid1, timvid1, "units", "years"
ncdf_attput, outid2, timvid2, "units", "years"
ncdf_attput, outid3, timvid3, "units", "years"
ncdf_attput, outid4, timvid4, "units", "years"
ncdf_attput, outid5, timvid5, "units", "years"
ncdf_attput, outid6, timvid6, "units", "years"
ncdf_attput, outid1, stfhtid, "units", "W/m^2"
ncdf_attput, outid2, stfsalid, "units", "s^-1"
ncdf_attput, outid3, tempid, "units", "degC"
ncdf_attput, outid4, salid, "units", "psu"
ncdf_attput, outid5, sstid, "units", "degC"
ncdf_attput, outid6, sssid, "units", "psu"
ncdf_attput, outid1, stfhtid, "missing_value", missing_value
ncdf_attput, outid2, stfsalid, "missing_value", missing_value
ncdf_attput, outid3, tempid, "missing_value", missing_value
ncdf_attput, outid4, salid, "missing_value", missing_value
ncdf_attput, outid5, sstid, "missing_value", missing_value
ncdf_attput, outid6, sssid, "missing_value", missing_value
ncdf_attput, outid1, lonvid1, "modulo", " "
ncdf_attput, outid2, lonvid2, "modulo", " "
ncdf_attput, outid3, lonvid3, "modulo", " "
ncdf_attput, outid4, lonvid4, "modulo", " "
ncdf_attput, outid5, lonvid5, "modulo", " "
ncdf_attput, outid6, lonvid6, "modulo", " "
ncdf_attput, outid1, lonvid1, "point_spacing", "even"
ncdf_attput, outid2, lonvid2, "point_spacing", "even"
ncdf_attput, outid3, lonvid3, "point_spacing", "even"
ncdf_attput, outid4, lonvid4, "point_spacing", "even"
ncdf_attput, outid5, lonvid5, "point_spacing", "even"
ncdf_attput, outid6, lonvid6, "point_spacing", "even"
ncdf_attput, outid1, latvid1, "point_spacing", "uneven"
ncdf_attput, outid2, latvid2, "point_spacing", "uneven"
ncdf_attput, outid3, latvid3, "point_spacing", "uneven"
ncdf_attput, outid4, latvid4, "point_spacing", "uneven"
ncdf_attput, outid5, latvid5, "point_spacing", "uneven"
ncdf_attput, outid6, latvid6, "point_spacing", "uneven"
ncdf_attput, outid3, depvid3, "point_spacing", "uneven"
ncdf_attput, outid4, depvid4, "point_spacing", "uneven"
ncdf_attput, outid1, timvid1, "point_spacing", "even"
ncdf_attput, outid2, timvid2, "point_spacing", "even"
ncdf_attput, outid3, timvid3, "point_spacing", "even"
ncdf_attput, outid4, timvid4, "point_spacing", "even"
ncdf_attput, outid5, timvid5, "point_spacing", "even"
ncdf_attput, outid6, timvid6, "point_spacing", "even"
ncdf_attput, outid3, depvid3, "positive", "down"
ncdf_attput, outid4, depvid4, "positive", "down"
ncdf_attput, outid1, /global, "title", title
ncdf_attput, outid2, /global, "title", title
ncdf_attput, outid3, /global, "title", title
ncdf_attput, outid4, /global, "title", title
ncdf_attput, outid5, /global, "title", title
ncdf_attput, outid6, /global, "title", title

; (5) Exit define mode
ncdf_control, outid1, /endef
ncdf_control, outid2, /endef
ncdf_control, outid3, /endef
ncdf_control, outid4, /endef
ncdf_control, outid5, /endef
ncdf_control, outid6, /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")
    depinid = ncdf_varid(inid, "zts")
    ncdf_varget, inid, loninid, lon
    ncdf_varget, inid, latinid, lat
    ncdf_varget, inid, depinid, z
    ncdf_varput, outid1, lonvid1, lon
    ncdf_varput, outid2, lonvid2, lon
    ncdf_varput, outid3, lonvid3, lon
    ncdf_varput, outid4, lonvid4, lon
    ncdf_varput, outid5, lonvid5, lon
    ncdf_varput, outid6, lonvid6, lon
    ncdf_varput, outid1, latvid1, lat
    ncdf_varput, outid2, latvid2, lat
    ncdf_varput, outid3, latvid3, lat
    ncdf_varput, outid4, latvid4, lat
    ncdf_varput, outid5, latvid5, lat
    ncdf_varput, outid6, latvid6, lat
    ncdf_varput, outid3, depvid3, z
    ncdf_varput, outid4, depvid4, z
    ncdf_varput, outid1, timvid1, yrs
    ncdf_varput, outid2, timvid2, yrs
    ncdf_varput, outid3, timvid3, yrs
    ncdf_varput, outid4, timvid4, yrs
    ncdf_varput, outid5, timvid5, yrs
    ncdf_varput, outid6, timvid6, yrs
  endif

  ; Read data from input file
  stfhtinid = ncdf_varid(inid, "stfht")
  stfsalinid = ncdf_varid(inid, "stfsal")
  tempinid = ncdf_varid(inid, "temp")
  salinid = ncdf_varid(inid, "sal")
  ncdf_varget, inid, stfhtinid, stfht
  ncdf_varget, inid, stfsalinid, stfsal
  ncdf_varget, inid, tempinid, temp
  ncdf_varget, inid, salinid, sal

  ; Write data to output file
  ncdf_varput, outid1, stfhtid, stfht, count=[nx, ny, 50], $
               offset=[0, 0, year-year1]
  ncdf_varput, outid2, stfsalid, stfsal, count=[nx, ny, 50], $
               offset=[0, 0, year-year1]
  ncdf_varput, outid3, tempid, temp, count=[nx, ny, nz, 50], $
               offset=[0, 0, 0, year-year1]
  ncdf_varput, outid4, salid, sal, count=[nx, ny, nz, 50], $
               offset=[0, 0, 0, year-year1]
  ncdf_varput, outid5, sstid, reform(temp(*, *, 0, *)), count=[nx, ny, 50], $
               offset=[0, 0, year-year1]
  ncdf_varput, outid6, sssid, reform(sal(*, *, 0, *)), count=[nx, ny, 50], $
               offset=[0, 0, year-year1]

endfor

; Close output files
ncdf_close, outid1
ncdf_close, outid2
ncdf_close, outid3
ncdf_close, outid4
ncdf_close, outid5
ncdf_close, outid6

end