pro enso_monthly

; Purpose
; -------
; Reads in monthly-mean MSLP and SST anomalies for a CSIRO model run, and
; derives the following quantities:
;
; (1) The Southern Oscillation Index (SOI)
; (2) The Nino 3.4 Index
; (3) The correlation between these two indices
; (4) The correlation between the SST at each gridpoint and the SOI
;
; If dP_T and dP_D are the MSLP anomalies at Tahiti (149deg34' W, 17deg32' S)
; and Darwin (130deg51' E, 12deg28' S) respectively, then the SOI is given by
;
;   SOI = 10.0 * (dp_T - dP_D) / SD(dp_T - dP_D)
;
; where SD() is the standard deviation. A separate standard deviation is
; calculated for each calendar month.
;
; The Nino 3.4 Index is simply the SST anomaly in the region 170-120 degW,
; 5 degS - 5 degN.
;
; History
; -------
; 2005 Jul 11	Steve Phipps	Original version, based on enso_annual.pro
; 2005 Aug 23	Steve Phipps	Modified to read in MSLP and SST anomalies,
;				rather than calculating them directly

; Define constants
nx = 64
ny = 56
deg_to_rad = !pi / 180.0
missing_value = -1.0e34

; Set longitudes and latitudes of Tahiti and Darwin
lont = 360.0 - (149.0 + 34.0/60.0)
latt = 0.0 - (17.0 + 32.0/60.0)
lond = 130.0 + 51.0/60.0
latd = 0.0 - (12.0 + 28.0/60.0)

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

; Get monthly-mean MSLP anomalies
filename = "spsl_" + run + "_mth_detrend.nc"
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, "psl_high")
lonid = ncdf_varid(ncid, "longitude")
latid = ncdf_varid(ncid, "latitude")
yrid = ncdf_varid(ncid, "year")
ncdf_varget, ncid, datid, psl_in
ncdf_varget, ncid, lonid, lon
ncdf_varget, ncid, latid, lat
ncdf_varget, ncid, yrid, years
ncdf_close, ncid

; Get monthly-mean SST anomalies
year1 = years(0)
year2 = years(n_elements(years)-1)
year1s = strtrim(string(year1), 1)
year2s = strtrim(string(year2), 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 + ".mth.detrend." + year1s + "-" + year2s + ".nc"
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, "sst_high")
lonid = ncdf_varid(ncid, "longitude")
latid = ncdf_varid(ncid, "latitude")
ncdf_varget, ncid, datid, sst_in
ncdf_varget, ncid, lonid, lono
ncdf_varget, ncid, latid, lato
ncdf_close, ncid

; Get other info from user
year1 = 0
year2 = 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, ""

; Declare output arrays
num_years = year2 - year1 +1
num_months = 12 * num_years
psl = fltarr(nx, ny, num_months)
sst = fltarr(nx, ny, num_months)
soi = fltarr(num_months)
soi5 = fltarr(num_months)
corr_sst_soi = fltarr(nx, ny)
nino34 = fltarr(num_months)
nino345 = fltarr(num_months)

; Move the MSLP and SST anomalies to 3-D arrays
for year = 0, num_years-1 do begin
  iyear = year + year1 - years(0)
  for month = 0, 11 do begin
    k = 12 * year + month
    psl(*, *, k) = psl_in(*, *, month, iyear)
    sst(*, *, k) = sst_in(*, *, month, iyear)
  endfor
endfor

; Derive monthly-mean MSLP anomalies at Tahiti
for i = 0, nx-2 do begin
  if ((lont gt lon(i)) and (lont le lon(i+1))) then begin
    imin = i
    weightx = 1.0 - (lont - lon(i)) / (lon(i+1) - lon(i))
  endif
endfor
for j = 0, ny-2 do begin
  if ((latt gt lat(j)) and (latt le lat(j+1))) then begin
    jmin = j
    weighty = 1.0 - (latt - lat(j)) / (lat(j+1) - lat(j))
  endif
endfor
dp_tahiti = weightx * weighty * psl(imin, jmin, *) + $
            (1.0-weightx) * weighty * psl(imin+1, jmin, *) + $
            weightx * (1.0-weighty) * psl(imin, jmin+1, *) + $
            (1.0-weightx) * (1.0-weighty) * psl(imin+1, jmin+1, *)
dp_tahiti = reform(dp_tahiti)

; Derive monthly-mean MSLP anomalies at Darwin
for i = 0, nx-2 do begin
  if ((lond gt lon(i)) and (lond le lon(i+1))) then begin
    imin = i
    weightx = 1.0 - (lond - lon(i)) / (lon(i+1) - lon(i))
  endif
endfor
for j = 0, ny-2 do begin
  if ((latd gt lat(j)) and (latd le lat(j+1))) then begin
    jmin = j
    weighty = 1.0 - (latd - lat(j)) / (lat(j+1) - lat(j))
  endif
endfor
dp_darwin = weightx * weighty * psl(imin, jmin, *) + $
            (1.0-weightx) * weighty * psl(imin+1, jmin, *) + $
            weightx * (1.0-weighty) * psl(imin, jmin+1, *) + $
            (1.0-weightx) * (1.0-weighty) * psl(imin+1, jmin+1, *)
dp_darwin = reform(dp_darwin)

; Derive SOI
dp = dp_tahiti - dp_darwin
for month = 0, 11 do begin
  index = 12 * indgen(num_years) + month
  soi(index) = 10.0 * dp(index) / stddev(dp(index))
endfor

; Derive Nino 3.4 Index. This bit of interpolation code is a bit lazy: it only
; works because the Nino 3.4 region doesn't span the Greenwich Meridian.
nino34(*) = 0.0
total_weight = 0.0
for j = 1, ny-2 do begin
  lats = 0.5 * (lato(j-1) + lato(j))
  latn = 0.5 * (lato(j) + lato(j+1))
  if ((latn lt -5.0) or (lats gt 5.0)) then continue
  weighty = 1.0
  if (lats lt -5.0) then weighty = $
    (sin(deg_to_rad * latn) + sin(deg_to_rad * 5.0)) / $
    (sin(deg_to_rad * latn) - sin(deg_to_rad * lats))
  if (latn gt 5.0) then weighty = $
    (sin(deg_to_rad * 5.0) - sin(deg_to_rad * lats)) / $
    (sin(deg_to_rad * latn) - sin(deg_to_rad * lats))
  for i = 1, nx-2 do begin
    lonw = 0.5 * (lono(i-1) + lono(i))
    lone = 0.5 * (lono(i) + lono(i+1))
    if ((lone lt 190.0) or (lonw gt 240.0)) then continue
    weightx = 1.0
    if (lonw lt 190.0) then weightx = (lone - 190.0) / (lone - lonw)
    if (lone gt 240.0) then weightx = (240.0 - lonw) / (lone - lonw)
    nino34 = nino34 + weightx * weighty * sst(i, j, *)
    total_weight = total_weight + weightx * weighty
  endfor
endfor
nino34 = nino34 / total_weight

; Smooth SOI and Nino 3.4 Index, by taking the five-month running mean
soi5(*) = missing_value
nino345(*) = missing_value
for i = 2, num_months-3 do begin
  soi5(i) = total(soi(i-2:i+2)) / 5.0
  nino345(i) = total(nino34(i-2:i+2)) / 5.0
endfor

; Derive correlations between SST and SOI
for j = 0, ny-1 do begin
  for i = 0, nx-1 do begin
    if (sst(i, j, 0, 0) lt -9000.0) then begin
      corr_sst_soi(i, j) = missing_value
    endif else begin
      corr_sst_soi(i, j) = correlate(sst(i, j, 2:num_months-3), $
                                     soi5(2:num_months-3))
    endelse
  endfor
endfor

; Derive correlation between Nino 3.4 index and SOI
corr_nino34_soi = correlate(nino345(2:num_months-3), soi5(2:num_months-3))
print, ""
print, "Correlation between Nino 3.4 Index and SOI = ", corr_nino34_soi
print, ""

; GENERATE OUTPUT FILE
;
; (1) Create netCDF file
year1s = strtrim(string(year1), 1)
year2s = strtrim(string(year2), 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 = "enso_mth_" + run + "_" + year1s + "-" + year2s + ".nc"
ncid = ncdf_create(filename)

; (2) Define dimensions
londid = ncdf_dimdef(ncid, "longitude", nx)
latdid = ncdf_dimdef(ncid, "latitude", ny)
mthdid = ncdf_dimdef(ncid, "month", num_months)

; (3) Define variables
lonvid = ncdf_vardef(ncid, "longitude", londid, /float)
latvid = ncdf_vardef(ncid, "latitude", latdid, /float)
mthvid = ncdf_vardef(ncid, "month", mthdid, /long)
dp_tahiti_id = ncdf_vardef(ncid, "dp_tahiti", mthdid, /float)
dp_darwin_id = ncdf_vardef(ncid, "dp_darwin", mthdid, /float)
soi_id = ncdf_vardef(ncid, "soi", mthdid, /float)
nino34_id = ncdf_vardef(ncid, "nino34", mthdid, /float)
corr_sst_soi_id = ncdf_vardef(ncid, "corr_sst_soi", [londid, latdid], /float)
corr_nino34_soi_id = ncdf_vardef(ncid, "corr_nino34_soi", /float)

; (4) Put attributes
ncdf_attput, ncid, dp_tahiti_id, "long_name", "MSLP anomaly at Tahiti"
ncdf_attput, ncid, dp_darwin_id, "long_name", "MSLP anomaly at Darwin"
ncdf_attput, ncid, soi_id, "long_name", "Southern Oscillation Index"
ncdf_attput, ncid, nino34_id, "long_name", "Nino 3.4 Index"
ncdf_attput, ncid, corr_sst_soi_id, "long_name", $
  "Correlation between SST and SOI"
ncdf_attput, ncid, corr_nino34_soi_id, "long_name", $
  "Correlation between Nino 3.4 Index and SOI"
ncdf_attput, ncid, lonvid, "units", "degrees_east"
ncdf_attput, ncid, latvid, "units", "degrees_north"
ncdf_attput, ncid, mthvid, "units", "months"
ncdf_attput, ncid, dp_tahiti_id, "units", "hPa"
ncdf_attput, ncid, dp_darwin_id, "units", "hPa"
ncdf_attput, ncid, nino34_id, "units", "K"
ncdf_attput, ncid, soi_id, "missing_value", missing_value
ncdf_attput, ncid, nino34_id, "missing_value", missing_value
ncdf_attput, ncid, corr_sst_soi_id, "missing_value", missing_value
ncdf_attput, ncid, lonvid, "point_spacing", "even"
ncdf_attput, ncid, latvid, "point_spacing", "uneven"
ncdf_attput, ncid, mthvid, "point_spacing", "even"
ncdf_attput, ncid, dp_tahiti_id, "longitude", lont
ncdf_attput, ncid, dp_tahiti_id, "latitude", latt
ncdf_attput, ncid, dp_darwin_id, "longitude", lond
ncdf_attput, ncid, dp_darwin_id, "latitude", latd
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, mthvid, 12*year1 + indgen(num_months) - 11
ncdf_varput, ncid, dp_tahiti_id, dp_tahiti
ncdf_varput, ncid, dp_darwin_id, dp_darwin
ncdf_varput, ncid, soi_id, soi5
ncdf_varput, ncid, nino34_id, nino345
ncdf_varput, ncid, corr_sst_soi_id, corr_sst_soi
ncdf_varput, ncid, corr_nino34_soi_id, corr_nino34_soi

; (7) Close file
ncdf_close, ncid

end