pro tsc_pcs_detrend

; Purpose
; -------
; Derives the principal components of annual-mean SAT 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 SAT anomalies that have been derived through the
; application of a high-pass filter.
;
; History
; -------
; 2005 Aug 24	Steve Phipps	Original version, based on tsc_pcs_area.pro

; Define constants
nx = 64
npcs = 10   ; Number of PCs to save to output file

; Get run name from user
run = ""
print, ""
read, run, prompt="Run name [e.g. d60]  "
print, ""

; Get annual-mean SAT anomalies
filename = "stsc_" + run + "_detrend.nc"
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, "tsc_high")
lonid = ncdf_varid(ncid, "longitude")
latid = ncdf_varid(ncid, "latitude")
yrid = ncdf_varid(ncid, "year")
ncdf_varget, ncid, datid, tsc_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_agcm_area.nc")
datid = ncdf_varid(ncid, "area")
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
tsc = double(tsc_in(*, jmin:jmax, year1-years(0):year2-years(0)))
lat = lat_in(jmin:jmax)
area = area_in(*, jmin:jmax)

; Declare output arrays
tsc_2d = dblarr(nx*ny, num_years)
tsc_pcs = dblarr(nx, ny, nx*ny)

; Normalise the data at each gridpoint, by subtracting the long-term mean
tsc_ltm = total(tsc, 3) / double(num_years)
for year = 0, num_years-1 do begin
  tsc(*, *, year) = tsc(*, *, year) - tsc_ltm
endfor

; Move the anomalies to a 2-D array, multiplying by the area of each gridbox
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    k = j * nx + i
    tsc_2d(k, *) = tsc(i, j, *) * area(i, j)
  endfor
endfor

; Generate the principal components
print, "Generating principal components ..."
tsc_trans = pcomp(tsc_2d, coefficients=tsc_pcs_2d, eigenvalues=eigenvalue, $
                  variances=variance, /covariance)

; Move the PCs back to a 3-D array, dividing by the area of each gridbox
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    k = j * nx + i
    tsc_pcs(i, j, *) = tsc_pcs_2d(k, *) / area(i, j)
  endfor
endfor

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

; Display some statistics
print, ""
print, "nx            = ", nx
print, "ny            = ", ny
print, "Number of PCs = ", nx*ny
print, ""
print, "Area analysed         = ", total(area), " m^2"
print, ""
print, "Sum of input data = ", total(tsc_2d), " K m^2"
print, ""
print, "Variance of input data       = ", total(tsc_2d^2) / $
                                          double(num_years-1), " K^2 m^4"
print, "Variance of transformed data = ", total(tsc_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 = "stsc_" + 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)
tsc_pcs_id = ncdf_vardef(ncid, "tsc_pcs", [londid, latdid, pcdid], /float)
tsc_trans_id = ncdf_vardef(ncid, "tsc_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, tsc_pcs_id, "long_name", "Principal components of SAT"
ncdf_attput, ncid, tsc_trans_id, "long_name", "Transformed SAT"
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, tsc_pcs_id, "units", "K m^2"
ncdf_attput, ncid, tsc_trans_id, "units", "K m^2"
ncdf_attput, ncid, eigenvalue_id, "units", "K^2 m^4"
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, tsc_pcs_id, tsc_pcs(*, *, 0:npcs-1)
ncdf_varput, ncid, tsc_trans_id, tsc_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