pro csiro_detrend_mth_climat

; Purpose
; -------
; Reads monthly-mean CSIRO AGCM output and generates two timeseries: a low-pass
; filtered timeseries, in which variability on timescales shorter than 100
; years has been removed, and a timeseries of the anomalies.
;
; History
; -------
; 2005 Aug 24	Steve Phipps	Original version, based on
;				csiro_detrend_climat.pro

; Define parameters
nx = 64
ny = 56

; Get various information from user
var = ""
run_name = ""
read, var, prompt="Enter the variable name [e.g. tsc]  "
read, run_name, prompt="Enter the run name      [e.g. d60]  "
print, ""

; Get monthly-mean model output
filename = "s" + var + "_" + run_name + ".nc"
ncid = ncdf_open(filename)
lonid = ncdf_varid(ncid, "longitude")
latid = ncdf_varid(ncid, "latitude")
yrid = ncdf_varid(ncid, "year")
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, lonid, lon
ncdf_varget, ncid, latid, lat
ncdf_varget, ncid, yrid, year
ncdf_varget, ncid, datid, data
ncdf_attget, ncid, datid, "long_name", long_name_in
long_name = string(long_name_in)
ncdf_attget, ncid, /global, "title", title_in
title = string(title_in)
num_years = n_elements(year)
print, "Variable = ", long_name
print, ""
print, "Data covers years ", strtrim(string(year(0)), 2), " to ", $
       strtrim(string(year(num_years-1)), 2), "."
print, ""
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
print, ""
ncdf_close, ncid

; Get other info from user
year1 = 0
year2 = 0
read, year1, prompt="First year to analyse   [e.g. 1]  "
read, year2, prompt="Last year to analyse [e.g. 1000]  "
print, ""

; Declare output arrays
l1 = year1-year(0)
l2 = year2-year(0)
num_years = year2 - year1 + 1
data_low = fltarr(nx, ny, 12, num_years)

; Apply a low-pass filter to derive the two timeseries
for month = 0, 11 do begin
  data_low(*, *, month, *) = lowpass_3d(reform(data(*, *, month, l1:l2)), 100)
endfor
data_high = data(*, *, *, l1:l2) - data_low

; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
year1s = strtrim(string(year1), 1)
year2s = strtrim(string(year2), 1)
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_name + "_mth_detrend_" + year1s + "-" + $
           year2s + ".nc"
ncid = ncdf_create(filename)

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

; (3) Define variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
mthvid = ncdf_vardef(ncid, "month", mthdid, /long)
yrvid = ncdf_vardef(ncid, "year", yrdid, /long)
lowid = ncdf_vardef(ncid, var+"_low", [londid, latdid, mthdid, yrdid], /float)
highid = ncdf_vardef(ncid, var+"_high", [londid, latdid, mthdid, yrdid], $
                     /float)

; (4) Put attributes
ncdf_attput, ncid, lowid, "long_name", "Low-pass filtered " + long_name
ncdf_attput, ncid, highid, "long_name", "Detrended " + long_name
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, mthvid, "units", "months"
ncdf_attput, ncid, yrvid, "units", "years"
if (has_units eq 1) then begin
  ncdf_attput, ncid, lowid, "units", units
  ncdf_attput, ncid, highid, "units", units
endif
ncdf_attput, ncid, lonvid, "modulo", " "
ncdf_attput, ncid, mthvid, "modulo", " "
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, mthvid, "point_spacing", "even"
ncdf_attput, ncid, yrvid, "point_spacing", "even"
ncdf_attput, ncid, /global, "title", title

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

; (6) Write data
ncdf_varput, ncid, lonvid, lon
ncdf_varput, ncid, latvid, lat
ncdf_varput, ncid, mthdid, indgen(12)+1
ncdf_varput, ncid, yrvid, year(l1:l2)
ncdf_varput, ncid, lowid, data_low
ncdf_varput, ncid, highid, data_high

; (7) Close file
ncdf_close, ncid

end