pro find_elnino

; Purpose
; -------
; Reads in values of the Nino 3.4 Index, and detects El Nino events.
;
; We use the definition of Trenberth (1975) i.e. that an El Nino event can be
; said to occur if the five-month running mean of the SST anomaly in the Nino
; 3.4 region [170-120degW, 5degS-5degN] exceeds 0.4K for at least six
; consecutive months. 
;
; History
; -------
; 2005 Jul 11	Steve Phipps	Original version
; 2005 Aug 30	Steve Phipps	Modified to calculate standard errors in the
;				means

; Define constants
max_events = 1000
months = ["Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", $
          "Oct", "Nov", "Dec"]

; Get info from user
run = ""
filename = ""
print, ""
read, run, prompt="Run name [e.g. d60]  "
print, ""
read, filename, prompt="Name of file containing Nino 3.4 index :  "
print, ""

; Get values of Nino 3.4 index
ncid = ncdf_open(filename)
datid = ncdf_varid(ncid, "nino34")
mthid = ncdf_varid(ncid, "month")
ncdf_varget, ncid, datid, nino34
ncdf_varget, ncid, mthid, month
ncdf_close, ncid

; Open output file
filename = "elnino_" + run + ".txt"
openw, 1, filename
printf, 1, "Event     Start date    End date   Length   Max anomaly"
printf, 1, "-----     ----------    --------   ------   -----------"

; Search for El Nino events
num_months = n_elements(month)
length = fltarr(max_events)
anomaly = fltarr(max_events)
num_events = 0
max_length = 0
max_anomaly = 0.0
elnino = 0
for i = 3, num_months-3 do begin   ; No Nino 3.4 Index for first/last 2 months

  if (elnino eq 0) then begin

    ; The last month was not part of an El Nino event, so check if this month
    ; is the beginning of an event
    if ((nino34(i-1) lt 0.4) and (nino34(i) ge 0.4)) then begin
      elnino = 1
      start_month = i
    endif

  endif else begin

    ; The last month was part of a possible El Nino event, so check if this
    ; month is the end of an event

    if ((nino34(i-1) ge 0.4) and (nino34(i) lt 0.4)) then begin

      ; The event is over, but it's only an El Nino event if it lasted for at
      ; least six months
      elnino = 0
      end_month = i - 1
      duration = end_month - start_month + 1
      if (duration ge 6) then begin

        ; It's an El Nino event!
        num_events = num_events + 1
        anom = max(nino34(start_month:end_month))
        length(num_events-1) = float(duration)
        anomaly(num_events-1) = anom
        if (duration gt max_length) then max_length = duration
        if (anom gt max_anomaly) then max_anomaly = anom
        year1 = (month(start_month) - 1) / 12 + 1
        year2 = (month(end_month) - 1) / 12 + 1
        month1 = (month(start_month) - 1) mod 12
        month2 = (month(end_month) - 1) mod 12
        year1s = strtrim(string(year1), 2)
        year2s = strtrim(string(year2), 2)
        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
        month1s = months(month1)
        month2s = months(month2)
        printf, 1, format='(i5, a10, a4, a10, a4, i7, f13.6)', num_events, $
                   year1s, month1s, year2s, month2s, duration, anom

      endif

    endif

  endelse

endfor

; Write summary to output file
tau_mean = float(num_months-4) / (12.0 * float(num_events))
tau_sd = tau_mean / sqrt(float(num_events))
l_mean = total(length(0:num_events-1)) / float(num_events)
l_sd = stddev(length(0:num_events-1)) / sqrt(float(num_events))
mag_mean = total(anomaly(0:num_events-1)) / float(num_events)
mag_sd = stddev(anomaly(0:num_events-1)) / sqrt(float(num_events))
printf, 1, ""
printf, 1, ""
printf, 1, "Average period = ", tau_mean, " +/- ", tau_sd, " years"
printf, 1, ""
printf, 1, "Average length = ", l_mean, " +/- ", l_sd, " months"
printf, 1, "Maximum length = ", max_length, " months"
printf, 1, ""
printf, 1, "Average magnitude = ", mag_mean, " +/- ", mag_sd, " K"
printf, 1, "Maximum magnitude = ", max_anomaly, " K"

; Close output file
close, 1

end