pro generate_dsss

; Purpose
; -------
; Reads in climatological SSSs and OGCM climatological SSSs, and calculates the
; corresponding SSS corrections required by the coupled model. Both annual-mean
; values, and effective values for the start of each month, are calculated.
;
; Note that the model interpolates between the values at the start and end of
; each month, in order to derive daily values. We apply a small correction to
; the estimated values for the start of each month, in order to ensure
; conservation of the annual mean.
;
; History
; -------
; 2004 Nov 23	Steve Phipps	Original version

; 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
dsann = fltarr(nx, ny)
dsann2 = fltarr(nx, ny)
ds_start = fltarr(nx, ny, 12)

; Get climatological SSSs
filename = ''
var = ''
print, ""
read, filename, prompt="Enter name of file containing climatological SSSs :  "
read, var, prompt="Enter name of variable containing SSS data :  "
print, ""
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, sss_clim
ncdf_close, ncid

; Get climatological SSSs from the OGCM
filename = ''
var = ''
read, filename, prompt="Enter name of file containing OGCM SSSs :  "
read, var, prompt="Enter name of variable containing SSS data :  "
print, ""
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, sss_ogcm
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 the monthly-mean SSS corrections
ds = sss_clim - sss_ogcm

; Fix missing values
ds(where(sss_ogcm lt -9000.0)) = 0.0

; Calculate annual-mean SSS corrections
dsann(*, *) = 0.0
for month = 0, 11 do begin
  dsann = dsann + ds(*, *, month) * days(month) / 365.0
endfor

; Average consecutive monthly values in order to obtain best estimates for
; the values at the start of each month
for month = 0, 11 do begin
  if (month eq 0) then begin
    ds_start(*, *, 0) = 0.5 * (ds(*, *, 11) + ds(*, *, 0))
  endif else begin
    ds_start(*, *, month) = 0.5 * (ds(*, *, month-1) + ds(*, *, month))
  endelse
endfor

; Calculate the error in the annual mean that arises when we interpolate
; between these values in order to obtain daily values
dsann2(*, *, *) = 0.0
for month = 0, 11 do begin
  mp1 = month + 1
  if (mp1 eq 12) then mp1 = 0
  dsann2 = dsann2 + 0.5 * (ds_start(*, *, month) + $
                           ds_start(*, *, mp1)) * days(month) / 365.0
endfor
dserr = dsann2 - dsann

; Correct the estimated values, by subtracting this error
for month = 0, 11 do begin
  ds_start(*, *, month) = ds_start(*, *, month) - dserr
endfor

; Save the SSS corrections to a netCDF file
;
; (1) Create netCDF file
filename = ''
read, filename, prompt="Enter name of SSS correction 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)
dsid = ncdf_vardef(ncid, "ds", [londid, latdid, mondid], /float)
ds2id = ncdf_vardef(ncid, "ds2", [londid, latdid, mondid], /float)
dsannid = ncdf_vardef(ncid, "dsann", [londid, latdid], /float)
dserrid = ncdf_vardef(ncid, "dserr", [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, dsid, ds
ncdf_varput, ncid, ds2id, ds_start
ncdf_varput, ncid, dsannid, dsann
ncdf_varput, ncid, dserrid, dserr

; (7) Close netCDF file
ncdf_close, ncid

end