pro generate_flux_adjustments

; Purpose
; -------
; Reads in AGCM and OGCM surface fluxes, derives the flux adjustments and saves
; them to a netCDF file. Both annual-mean and adjusted monthly values are
; calculated.
;
; The adjustments to the monthly values are necessary if we are to treat them
; as representing values for the midpoint of each month, and then interpolate
; between them in order to obtain daily values. This approach introduces a
; slight conservation error, in that the differing lengths of the months cause
; a slight failure to conserve the annual mean. An annual-mean correction is
; applied to the values at each gridpoint, in order to correct for this error.
;
; NOTE that the TS grid is used in all cases, although we should use the UV
; grid in the case of the surface stresses.
;
; History
; -------
; 2004 Mar 7	Steve Phipps	Original version
; 2004 Nov 19	Steve Phipps	Modified to: 
;				(1) Calculate the annual means, and to apply
;				    corrections to the monthly values.
;				(2) To save the data to a netCDF file, rather
;				    than directly to a binary auxiliary file.

; Define constants
nx = 64
ny = 56
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]

; Initialise variables
dfann = fltarr(nx, ny)
dfann2 = fltarr(nx, ny)
df2 = fltarr(nx, ny, 12)

; Get AGCM surface fluxes
filename = ''
var = ''
read, filename, prompt="Enter name of file containing AGCM surface fluxes: "
read, var, prompt="Enter name of variable containing flux data: "
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, fa
ncdf_close, ncid

; Get OGCM surface fluxes
filename = ''
var = ''
read, filename, prompt="Enter name of file containing OGCM surface fluxes: "
read, var, prompt="Enter name of variable containing flux data: "
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, fo
ncdf_close, ncid

; Get CSIRO OGCM grid
ncid = ncdf_open("/home/sjphipps/phd/csiro/data/csiro_ogcm_ts_landsea.nc")
datid = ncdf_varid(ncid, "lonts")
ncdf_varget, ncid, datid, lonts
datid = ncdf_varid(ncid, "latts")
ncdf_varget, ncid, datid, latts
ncdf_close, ncid

; Derive monthly-mean flux adjustments
df = fa - fo

; Calculate annual-mean flux adjustments
dfann(*, *) = 0.0
for month = 0, 11 do begin
  dfann = dfann + df(*, *, month) * days(month) / 365.0
endfor

; Calculate the annual-mean error that arises when we interpolate between
; the monthly means in order to obtain daily values
dfann2(*, *) = 0.0
for month = 0, 11 do begin
  monthm1 = month - 1
  monthp1 = month + 1
  if (monthm1 eq -1) then monthm1 = 11
  if (monthp1 eq 12) then monthp1 = 0
  dfann2 = dfann2 + df(*, *, month) * (0.125 * days(monthm1) + $
                                       0.75 * days(month) + $
                                       0.125 * days(monthp1)) / 365.0
endfor
err = dfann2 - dfann

; Correct the monthly means, by subtracting this error
for month = 0, 11 do begin
  df2(*, *, month) = df(*, *, month) - err
endfor

; Save the flux adjustments to a netCDF file
;
; (1) Create netCDF file
filename = ''
read, filename, prompt="Enter name of flux adjustment file: "
ncid = ncdf_create(filename)

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

; (3) Declare variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
monvid = ncdf_vardef(ncid, "month", mondid, /long)
dfid = ncdf_vardef(ncid, "df", [londid, latdid, mondid], /float)
df2id = ncdf_vardef(ncid, "df2", [londid, latdid, mondid], /float)
dfannid = ncdf_vardef(ncid, "dfann", [londid, latdid], /float)
errid = ncdf_vardef(ncid, "error", [londid, latdid], /float)

; (4) Put attributes
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, monvid, "units", "months"
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"

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

; (6) Write data
ncdf_varput, ncid, lonvid, lonts
ncdf_varput, ncid, latvid, latts
ncdf_varput, ncid, monvid, indgen(12) + 1
ncdf_varput, ncid, dfid, df
ncdf_varput, ncid, df2id, df2
ncdf_varput, ncid, dfannid, dfann
ncdf_varput, ncid, errid, err

; (7) Close netCDF file
ncdf_close, ncid

end