pro csiro_climat_agcm_stats_multi

; Purpose
; -------
; Reads in monthly average values from CSIRO AGCM output and writes various
; climatological values, including the standard errors in each of the means, to
; a new netCDF file.
;
; History
; -------
; 2005 Nov 26	Steven Phipps	Original version, based on
;				csiro_climat_agcm_stats.pro and
;				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
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
var = ""
run_name = ""
year1 = 0
year2 = 0
read, var, prompt="Please enter the variable name [e.g. tsc]  "
read, run_name, prompt="Please enter the run name      [e.g. d73]  "
read, year1, prompt="Enter the initial year        [e.g. 1301]  "
read, year2, prompt="Enter the final year          [e.g. 1400]  "

; 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 array to hold data
data = fltarr(nx, ny, 12, num_years)

; 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 = "s" + var + "_" + run_name + "_" + y1s + "-" + y2s + ".nc"
  print, "Reading data from file ", filename, " ..."
  ncid = ncdf_open(filename)

  ; Read data
  datid = ncdf_varid(ncid, var)
  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
    ncdf_attget, ncid, datid, "long_name", long_name_in
    long_name = string(long_name_in)
    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
  endif

  ; Close netCDF file
  ncdf_close, ncid

  ; Fix missing values
  index = where(data gt 9.0e36)
  if (index ne -1) then data_in(index) = missing_value

  ; Add data to main array
  data(*, *, *, y1-year1:y2-year1) = data_in

endfor

; Derive climatology
clim_month = fltarr(nx, ny, 12, 2)
clim_ann = fltarr(nx, ny, 2)
clim_djf = fltarr(nx, ny, 2)
clim_jja = fltarr(nx, ny, 2)
ann_temp = fltarr(num_years)
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin

    ; Monthly climatology
    for k = 0, 11 do begin
      temp = data(i, j, k, *)
      if (min(temp) lt -1.0e30) then begin
        clim_month(i, j, k, *) = missing_value
      endif else begin
        stats = moment(temp)
        clim_month(i, j, k, 0) = stats(0)
        clim_month(i, j, k, 1) = sqrt(stats(1) / float(num_years))
      endelse
    endfor

    ; Annual means
    ann_temp(*) = 0.0
    for k = 0, 11 do begin
      ann_temp = ann_temp + days(k) * data(i, j, k, *) / 365.0
    endfor
    if (min(ann_temp) lt -1.0e30) then begin
      clim_ann(i, j, *) = missing_value
    endif else begin
      stats = moment(ann_temp)
      clim_ann(i, j, 0) = stats(0)
      clim_ann(i, j, 1) = sqrt(stats(1) / float(num_years))
    endelse

    ; DJF means
    djf_temp = (days(11)*data(i, j, 11, *) + $
                days(0)*data(i, j, 0, *) + $
                days(1)*data(i, j, 1, *)) / 90.0
    if (min(djf_temp) lt -1.0e30) then begin
      clim_djf(i, j, *) = missing_value
    endif else begin
      stats = moment(djf_temp)
      clim_djf(i, j, 0) = stats(0)
      clim_djf(i, j, 1) = sqrt(stats(1) / float(num_years))
    endelse

    ; JJA means
    jja_temp = (days(5)*data(i, j, 5, *) + $
                days(6)*data(i, j, 6, *) + $
                days(7)*data(i, j, 7, *)) / 92.0
    if (min(jja_temp) lt -1.0e30) then begin
      clim_jja(i, j, *) = missing_value
    endif else begin
      stats = moment(jja_temp)
      clim_jja(i, j, 0) = stats(0)
      clim_jja(i, j, 1) = sqrt(stats(1) / float(num_years))
    endelse

  endfor
endfor

; Derive name of output 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 + "_climats_" + year1s + "-" + year2s + $
           ".nc"
print, "Output file = ", filename

; Create new netCDF file
ncid = ncdf_create(filename)

; Declare dimensions
londid = ncdf_dimdef(ncid, "longitude", nx)
latdid = ncdf_dimdef(ncid, "latitude", ny)
mondid = ncdf_dimdef(ncid, "month", 12)

; Declare variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
monvid = ncdf_vardef(ncid, "month", mondid, /long)
climid = ncdf_vardef(ncid, var, [londid, latdid, mondid], /float)
annid = ncdf_vardef(ncid, var+"ann", [londid, latdid], /float)
djfid = ncdf_vardef(ncid, var+"djf", [londid, latdid], /float)
jjaid = ncdf_vardef(ncid, var+"jja", [londid, latdid], /float)
climseid = ncdf_vardef(ncid, var+"se", [londid, latdid, mondid], /float)
annseid = ncdf_vardef(ncid, var+"annse", [londid, latdid], /float)
djfseid = ncdf_vardef(ncid, var+"djfse", [londid, latdid], /float)
jjaseid = ncdf_vardef(ncid, var+"jjase", [londid, latdid], /float)

; Put attributes
ncdf_attput, ncid, climid, "long_name", long_name
ncdf_attput, ncid, annid, "long_name", "Annual-mean " + long_name
ncdf_attput, ncid, djfid, "long_name", "DJF " + long_name
ncdf_attput, ncid, jjaid, "long_name", "JJA " + long_name
ncdf_attput, ncid, climseid, "long_name", "Standard error of " + long_name
ncdf_attput, ncid, annseid, "long_name", "Standard error of annual-mean " + $
                                         long_name
ncdf_attput, ncid, djfseid, "long_name", "Standard error of DJF " + long_name
ncdf_attput, ncid, jjaseid, "long_name", "Standard error of JJA " + long_name
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, monvid, "units", "months"
if (has_units eq 1) then begin
  ncdf_attput, ncid, climid, "units", units
  ncdf_attput, ncid, annid, "units", units
  ncdf_attput, ncid, djfid, "units", units
  ncdf_attput, ncid, jjaid, "units", units
  ncdf_attput, ncid, climseid, "units", units
  ncdf_attput, ncid, annseid, "units", units
  ncdf_attput, ncid, djfseid, "units", units
  ncdf_attput, ncid, jjaseid, "units", units
endif
ncdf_attput, ncid, climid, "missing_value", missing_value
ncdf_attput, ncid, annid, "missing_value", missing_value
ncdf_attput, ncid, djfid, "missing_value", missing_value
ncdf_attput, ncid, jjaid, "missing_value", missing_value
ncdf_attput, ncid, climseid, "missing_value", missing_value
ncdf_attput, ncid, annseid, "missing_value", missing_value
ncdf_attput, ncid, djfseid, "missing_value", missing_value
ncdf_attput, ncid, jjaseid, "missing_value", missing_value
ncdf_attput, ncid, lonvid, "modulo", " "
ncdf_attput, ncid, monvid, "modulo", " "
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, monvid, "point_spacing", "even"
ncdf_attput, ncid, /global, "title", "CSIRO model run " + run_name + $
  ", years " + year1s + " to " + year2s

; Exit define mode
ncdf_control, ncid, /endef

; Write data
ncdf_varput, ncid, lonvid, lon
ncdf_varput, ncid, latvid, lat
ncdf_varput, ncid, monvid, indgen(12) + 1
ncdf_varput, ncid, climid, clim_month(*, *, *, 0)
ncdf_varput, ncid, annid, clim_ann(*, *, 0)
ncdf_varput, ncid, djfid, clim_djf(*, *, 0)
ncdf_varput, ncid, jjaid, clim_jja(*, *, 0)
ncdf_varput, ncid, climseid, clim_month(*, *, *, 1)
ncdf_varput, ncid, annseid, clim_ann(*, *, 1)
ncdf_varput, ncid, djfseid, clim_djf(*, *, 1)
ncdf_varput, ncid, jjaseid, clim_jja(*, *, 1)

; Close netCDF file
ncdf_close, ncid

end