pro average_qflux

; Purpose
; -------
; Reads in monthly qflux output from the CSIRO AGCM, derives a climatology
; and saves the values to file.
;
; History
; -------
; 2002 Nov 8	Steven Phipps	Original version
; 2003 May 24	Steven Phipps	Modified to reflect use of 32-bit integers by
;				model, rather than 64-bit as previously
; 2009 May 2	Steven Phipps	Modified for five-digit year numbers

; Define constants
lon = 64
lat = 28
nmonth = 12

; Get run details from user
run = ""
year1 = 0
year2 = 0
read, run, prompt="Run name [e.g. b44]  "
read, year1, prompt="Initial year for climatology [e.g. 6]   "
read, year2, prompt="Final year for climatology   [e.g. 35]  "

; Declare arrays to hold data
nyear = year2 - year1 + 1
ochf = dblarr(lon, lat, 2, nmonth, nyear)
occur = dblarr(lon, lat, 2, nmonth, nyear)
ochf_in = dblarr(lon, 2, lat)
occur_in = dblarr(lon, 2, lat)
exptyp = bytarr(50)
kday = 0l
kdayp = 0l
chtst = bytarr(4)

; Loop over years
for year = year1, year2 do begin
  years = strtrim(string(year), 2)
  if (year lt 10000) then years = '0' + years
  if (year lt 1000) then years = '0' + years
  if (year lt 100) then years = '0' + years
  if (year lt 10) then years = '0' + years

  ; Loop over months
  for month = 1, 12 do begin
    months = strtrim(string(month), 2)
    if (month lt 10) then months = '0' + months

    ; Open file for unformatted read access
    filename = 'qflx' + years + months + '.' + run
    openr, 1, filename, /f77_unformatted

    ; Read data
    readu, 1, exptyp
    readu, 1, kday, kdayp, chtst
    readu, 1, ochf_in
    readu, 1, occur_in

    ; Add data to output arrays
    ochf(*, *, 0, month-1, year-year1) = ochf_in(*, 0, *)
    ochf(*, *, 1, month-1, year-year1) = ochf_in(*, 1, *)
    occur(*, *, 0, month-1, year-year1) = occur_in(*, 0, *)
    occur(*, *, 1, month-1, year-year1) = occur_in(*, 1, *)

    ; Close file
    close, 1

  endfor
endfor

; Derive climatologies
ochf_clim = total(ochf, 5) / float(nyear)
occur_clim = total(occur, 5) / float(nyear)

; Save data
outfile1 = 'qflux.dat.' + run
outfile2 = 'qfluxform.dat.' + run
outfile3 = 'ochfform.dat.' + run
openw, 1, outfile1, /f77_unformatted
openw, 2, outfile2
openw, 3, outfile3
for month = 1, 12 do begin
  header = string(replicate(32b, 80))
  title = 'Q-flux: run ' + run + ', average for month ' + strtrim(string(month), 2)
  strput, header, title
  writeu, 1, header
  printf, 2, header
  printf, 3, header
  writeu, 1, long(month)
  printf, 2, long(month)
  printf, 3, long(month)
  writeu, 1, occur_clim(*, *, *, month-1)
  printf, 2, format='(6f15.8)', occur_clim(*, *, *, month-1)
  printf, 3, format='(6f15.8)', ochf_clim(*, *, *, month-1)
endfor
close, 1, 2, 3

end