pro sst_pcs_detrend

; Purpose
; -------
; Derives the principal components of annual-mean SST for a CSIRO model run.
;
; This version weights the values by the area of each gridbox, ensuring that
; physically-correct principal components are derived.
;
; This version also uses SST anomalies that have been derived through the
; application of a high-pass filter.
;
; History
; -------
; 2005 Aug 23	Steve Phipps	Original version, based on sst_pcs_area.pro

; Define constants
nx = 64
npcs = 10   ; Number of PCs to save to output file
missing_value = -1.0e34

; Get information from user
run = ""
filename = ""
print, ""
read, run, prompt="Run name                              [e.g. d60]  "
read, filename, prompt="Input file [e.g. sst.d60.detrend.00001-01000.nc]  "
print, ""

; Get SST anomalies
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, "sst_high")
lonid = ncdf_varid(ncid, "longitude")
latid = ncdf_varid(ncid, "latitude")
yrid = ncdf_varid(ncid, "year")
ncdf_varget, ncid, datid, sst_in
ncdf_varget, ncid, lonid, lon
ncdf_varget, ncid, latid, lat_in
ncdf_varget, ncid, yrid, years
ncdf_close, ncid

; Get areas of gridboxes
ncid = ncdf_open("csiro_ogcm_ts_area_volume.nc")
datid = ncdf_varid(ncid, "da")
ncdf_varget, ncid, datid, area_in
ncdf_close, ncid

; Get other info from user
year1 = 0
year2 = 0
ny = 0
print, "Data covers years ", years(0), " to ", years(n_elements(years)-1), "."
print, ""
read, year1, prompt="First year to analyse   [e.g. 1]  "
read, year2, prompt="Last year to analyse [e.g. 1000]  "
print, ""
read, ny, prompt="Number of latitude rows to process [2-56]  "
print, ""

; Extract data to analyse
jmin = 28 - ny/2
jmax = 27 + ny/2
num_years = year2 - year1 + 1
sst = double(sst_in(*, jmin:jmax, year1-years(0):year2-years(0)))
lat = lat_in(jmin:jmax)
area = area_in(jmin:jmax)

; Derive the number of ocean gridpoints, and the total surface area
nocean = 0
total_area = 0.0d
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    if (sst(i, j, 0) gt -9000.0d) then begin
      nocean = nocean + 1
      total_area = total_area + area(j)
    endif
  endfor
endfor

; Derive arrays containing the co-ordinates of each ocean gridpoint
iocean = intarr(nocean)
jocean = intarr(nocean)
k = 0
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    if (sst(i, j, 0) gt -9000.0d) then begin
      iocean(k) = i
      jocean(k) = j
      k = k + 1
    endif
  endfor
endfor
if (k ne nocean) then begin
  print, "ERROR: Wrong number of ocean gridpoints!"
  stop
endif

; Declare output arrays
sst_2d = dblarr(nocean, num_years)
sst_pcs = dblarr(nx, ny, nocean)

; Normalise the data at each gridpoint, by subtracting the long-term mean
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    if (sst(i, j, 0) gt -9000.0d) then begin
      sst(i, j, *) = sst(i, j, *) - total(sst(i, j, *)) / double(num_years)
    endif
  endfor
endfor

; Move the anomalies to a 2-D array, multiplying by the area of each gridbox
for k = 0, nocean-1 do begin
  sst_2d(k, *) = sst(iocean(k), jocean(k), *) * area(jocean(k))
endfor

; Generate the principal components
print, "Generating principal components ..."
sst_trans = pcomp(sst_2d, coefficients=sst_pcs_2d, eigenvalues=eigenvalue, $
                  variances=variance, /covariance)

; Move the PCs back to a 3-D array, dividing by the area of each gridbox
sst_pcs(*, *, *) = missing_value
for k = 0, nocean-1 do begin
  sst_pcs(iocean(k), jocean(k), *) = sst_pcs_2d(k, *) / area(jocean(k))
endfor

; Rescale the transformed data, so that the coefficients represent multiples of
; normalised eigenvectors (i.e. eigenvectors with unit magnitude)
for k = 0, nocean-1 do begin
  sst_trans(k, *) = sst_trans(k, *) / sqrt(eigenvalue(k))
endfor

; Display some statistics
print, ""
print, "nx            = ", nx
print, "ny            = ", ny
print, "Number of PCs = ", nocean
print, ""
print, "Total surface area = ", total_area, " m^2"
print, ""
print, "Sum of input data = ", total(sst_2d), " K m^2"
print, ""
print, "Variance of input data       = ", total(sst_2d^2) / $
                                          double(num_years-1), " K^2 m^4"
print, "Variance of transformed data = ", total(sst_trans^2) / $
                                          double(num_years-1), " K^2 m^4"
print, "Sum of eigenvalues           = ", total(eigenvalue), " K^2 m^4"
print, ""
print, "Percent variance of first 10 PCs"
print, "--------------------------------"
for i = 0, 9 do begin
  print, i+1, 100.0d*variance(i)
endfor
print, ""
 
; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
year1s = strtrim(string(year1), 1)
year2s = strtrim(string(year2), 1)
nys = strtrim(string(ny), 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 = "sst." + run + ".pcs." + nys + "." + year1s + "-" + year2s + ".nc"
ncid = ncdf_create(filename)

; (2) Define dimensions
londid = ncdf_dimdef(ncid, "longitude", nx)
latdid = ncdf_dimdef(ncid, "latitude", ny)
yrdid = ncdf_dimdef(ncid, "year", num_years)
pcdid = ncdf_dimdef(ncid, "pc", npcs)

; (3) Define variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
yrvid = ncdf_vardef(ncid, "year", yrdid, /long)
pcvid = ncdf_vardef(ncid, "pc", pcdid, /long)
sst_pcs_id = ncdf_vardef(ncid, "sst_pcs", [londid, latdid, pcdid], /float)
sst_trans_id = ncdf_vardef(ncid, "sst_trans", [pcdid, yrdid], /float)
eigenvalue_id = ncdf_vardef(ncid, "eigenvalue", pcdid, /float)
variance_id = ncdf_vardef(ncid, "variance", pcdid, /float)

; (4) Put attributes
ncdf_attput, ncid, sst_pcs_id, "long_name", "Principal components of SST"
ncdf_attput, ncid, sst_trans_id, "long_name", "Transformed SST"
ncdf_attput, ncid, eigenvalue_id, "long_name", "Eigenvalue"
ncdf_attput, ncid, variance_id, "long_name", "Fraction of total variance"
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, yrvid, "units", "years"
ncdf_attput, ncid, sst_pcs_id, "units", "K m^2"
ncdf_attput, ncid, sst_trans_id, "units", "K m^2"
ncdf_attput, ncid, eigenvalue_id, "units", "K^2 m^4"
ncdf_attput, ncid, sst_pcs_id, "missing_value", missing_value
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, yrvid, "point_spacing", "even"
ncdf_attput, ncid, pcvid, "point_spacing", "even"
ncdf_attput, ncid, /global, "title", "CSIRO model run " + run + ", years " + $
                                     year1s + " to " + year2s

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

; (6) Write data
ncdf_varput, ncid, lonvid, lon
ncdf_varput, ncid, latvid, lat
ncdf_varput, ncid, yrvid, years(year1-years(0):year2-years(0))
ncdf_varput, ncid, pcvid, indgen(npcs) + 1
ncdf_varput, ncid, sst_pcs_id, sst_pcs(*, *, 0:npcs-1)
ncdf_varput, ncid, sst_trans_id, sst_trans(0:npcs-1, *)
ncdf_varput, ncid, eigenvalue_id, reform(eigenvalue(0:npcs-1))
ncdf_varput, ncid, variance_id, reform(variance(0:npcs-1))

; (7) Close file
ncdf_close, ncid

end