pro generate_dsst

; Purpose
; -------
; Reads in climatological SSTs, climatological values for DTM from the AGCM and
; climatological SSTs from the OGCM, and calculates the corresponding SST
; corrections required by the coupled model. Both annual-mean values, and
; effective values for the start of each month, are calculated, for both the
; SST corrections and DTM.
;
; 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 Mar 8	Steve Phipps	Original version
; 2004 Nov 19	Steve Phipps	Modified to:
;				(1) Interpolate the OGCM SSTs onto the AGCM
;				    grid, as in the coupled model.
;                               (2) Calculate the annual means, and to apply
;                                   corrections to the monthly values.
;                               (3) 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
dt = fltarr(nx, ny, 12)
dtann = fltarr(nx, ny)
dtmann = fltarr(nx, ny)
dtann2 = fltarr(nx, ny)
dtmann2 = fltarr(nx, ny)
dt_start = fltarr(nx, ny, 12)
dtm_start = fltarr(nx, ny, 12)

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

; Get climatological values for DTM from the AGCM
filename = ''
var = ''
read, filename, prompt="Enter name of file containing AGCM DTM :  "
read, var, prompt="Enter name of variable containing DTM data :  "
print, ""
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, dtm
ncdf_close, ncid

; Get climatological SSTs from the OGCM
filename = ''
var = ''
read, filename, prompt="Enter name of file containing OGCM SSTs :  "
read, var, prompt="Enter name of variable containing SST data :  "
print, ""
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, var)
ncdf_varget, ncid, datid, sst_ogcm
ncdf_close, ncid

; Get CSIRO OGCM grid and land/sea mask
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
datid = ncdf_varid(ncid, "mask")
ncdf_varget, ncid, datid, masko
ncdf_close, ncid

; Get areas of gridpoints on CSIRO OGCM grid
ncid = ncdf_open("/home/sjphipps/phd/csiro/data/csiro_ogcm_ts_area.nc")
datid = ncdf_varid(ncid, "area")
ncdf_varget, ncid, datid, area
ncdf_close, ncid

; Get CSIRO AGCM land/sea mask
ncid = ncdf_open("/home/sjphipps/phd/csiro/data/csiro_landsea.nc")
datid = ncdf_varid(ncid, "mask")
ncdf_varget, ncid, datid, maska
ncdf_close, ncid

; Interpolate OGCM SSTs onto the AGCM grid, deriving the monthly-mean SST
; corrections in the process
;
; (1) Estimate SSTs for those gridpoints which are regarded as being ocean by
;     the AGCM, but land by the OGCM.
sst_ogcm(26, 15, *) = 0.5 * (sst_ogcm(25, 15, *) + sst_ogcm(27, 15, *))
sst_ogcm(7, 20, *) = 0.75 * sst_ogcm(7, 19, *) + 0.25 * sst_ogcm(7, 23, *)
sst_ogcm(7, 21, *) = 0.5 * (sst_ogcm(7, 19, *) + sst_ogcm(7, 23, *))
sst_ogcm(7, 22, *) = 0.25 * sst_ogcm(7, 19, *) + 0.75 * sst_ogcm(7, 23, *)
sst_ogcm(26, 24, *) = 0.5 * (sst_ogcm(25, 24, *) + sst_ogcm(27, 24, *))
sst_ogcm(2, 26, *) = sst_ogcm(1, 26, *)
sst_ogcm(20, 27, *) = (2.0 * sst_ogcm(20, 26, *) + sst_ogcm(20, 29, *)) / 3.0
sst_ogcm(20, 28, *) = (sst_ogcm(20, 26, *) + 2.0 * sst_ogcm(20, 29, *)) / 3.0
sst_ogcm(19, 28, *) = sst_ogcm(19, 29, *)
sst_ogcm(18, 29, *) = 0.5 * (sst_ogcm(17, 29, *) + sst_ogcm(19, 29, *))
sst_ogcm(18, 30, *) = 0.5 * (sst_ogcm(17, 30, *) + sst_ogcm(19, 30, *))
sst_ogcm(18, 31, *) = sst_ogcm(17, 31, *)
sst_ogcm(22, 38, *) = 0.5 * (sst_ogcm(22, 37, *) + sst_ogcm(22, 39, *))
sst_ogcm(63, 39, *) = 0.5 * (sst_ogcm(62, 39, *) + sst_ogcm(0, 39, *))
sst_ogcm(51, 47, *) = (2.0 * sst_ogcm(50, 47, *) + sst_ogcm(53, 47, *)) / 3.0
sst_ogcm(52, 47, *) = (sst_ogcm(50, 47, *) + 2.0 * sst_ogcm(53, 47, *)) / 3.0
sst_ogcm(34, 48, *) = 0.5 * (sst_ogcm(34, 47, *) + sst_ogcm(34, 49, *))
sst_ogcm(51, 51, *) = 0.5 * (sst_ogcm(50, 51, *) + sst_ogcm(52, 51, *))

; (2) Ensure conservation of global-mean SST, deriving the monthly-mean SST
;     corrections in the process.
for month = 0, 11 do begin
  atoto = 0.0
  atota = 0.0
  ttoto = 0.0
  ttota = 0.0
  for j = 0, ny-1 do begin
    for i = 0, nx-1 do begin
      if (masko(i, j) eq 0) then begin
        atoto = atoto + area(i, j)
        ttoto = ttoto + area(i, j) * sst_ogcm(i, j, month)
      endif
      if (maska(i, j) eq 0) then begin
        atota = atota + area(i, j)
        ttota = ttota + area(i, j) * sst_ogcm(i, j, month)
      endif
    endfor
  endfor
  meano = ttoto / atoto
  meana = ttota / atota
  corr = meano - meana
  for j = 0, ny-1 do begin
    for i = 0, nx-1 do begin
      if (maska(i, j) eq 0) then begin
        dt(i, j, month) = sst_clim(i, j, month) + dtm(i, j, month) - $
                          sst_ogcm(i, j, month) - corr
      endif else begin
        dt(i, j, month) = 0.0
      endelse
    endfor
  endfor
  print, "Month = ", month+1
  print, ""
  print, "Ocean area on OGCM grid = ", atoto, " m^2"
  print, "Ocean area on AGCM grid = ", atota, " m^2"
  print, ""
  print, "Mean SST on OGCM grid = ", meano, " degC"
  print, "Mean SST on AGCM grid = ", meana, " degC"
  print, ""
  print, "Correction to SSTs on AGCM grid = ", corr, " degC"
  print, ""
endfor

; Calculate annual-mean SST corrections and DTM
dtann(*, *) = 0.0
dtmann(*, *) = 0.0
for month = 0, 11 do begin
  dtann = dtann + dt(*, *, month) * days(month) / 365.0
  dtmann = dtmann + dtm(*, *, 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
    dt_start(*, *, 0) = 0.5 * (dt(*, *, 11) + dt(*, *, 0))
    dtm_start(*, *, 0) = 0.5 * (dtm(*, *, 11) + dtm(*, *, 0))
  endif else begin
    dt_start(*, *, month) = 0.5 * (dt(*, *, month-1) + dt(*, *, month))
    dtm_start(*, *, month) = 0.5 * (dtm(*, *, month-1) + dtm(*, *, month))
  endelse
endfor

; Calculate the error in the annual mean that arises when we interpolate
; between these values in order to obtain daily values
dtann2(*, *, *) = 0.0
dtmann2(*, *, *) = 0.0
for month = 0, 11 do begin
  mp1 = month + 1
  if (mp1 eq 12) then mp1 = 0
  dtann2 = dtann2 + 0.5 * (dt_start(*, *, month) + $
                           dt_start(*, *, mp1)) * days(month) / 365.0
  dtmann2 = dtmann2 + 0.5 * (dtm_start(*, *, month) + $
                             dtm_start(*, *, mp1)) * days(month) / 365.0
endfor
dterr = dtann2 - dtann
dtmerr = dtmann2 - dtmann

; Correct the estimated values, by subtracting this error
for month = 0, 11 do begin
  dt_start(*, *, month) = dt_start(*, *, month) - dterr
  dtm_start(*, *, month) = dtm_start(*, *, month) - dtmerr
endfor

; Save the SST corrections to a netCDF file
;
; (1) Create netCDF file
filename = ''
read, filename, prompt="Enter name of SST 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)
dtid = ncdf_vardef(ncid, "dt", [londid, latdid, mondid], /float)
dtmid = ncdf_vardef(ncid, "dtm", [londid, latdid, mondid], /float)
dt2id = ncdf_vardef(ncid, "dt2", [londid, latdid, mondid], /float)
dtm2id = ncdf_vardef(ncid, "dtm2", [londid, latdid, mondid], /float)
dtannid = ncdf_vardef(ncid, "dtann", [londid, latdid], /float)
dtmannid = ncdf_vardef(ncid, "dtmann", [londid, latdid], /float)
dterrid = ncdf_vardef(ncid, "dterr", [londid, latdid], /float)
dtmerrid = ncdf_vardef(ncid, "dtmerr", [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, dtid, dt
ncdf_varput, ncid, dtmid, dtm
ncdf_varput, ncid, dt2id, dt_start
ncdf_varput, ncid, dtm2id, dtm_start
ncdf_varput, ncid, dtannid, dtann
ncdf_varput, ncid, dtmannid, dtmann
ncdf_varput, ncid, dterrid, dterr
ncdf_varput, ncid, dtmerrid, dtmerr

; (7) Close netCDF file
ncdf_close, ncid

end