;======================================
; AinsworthEtAl_2016_gbrReefs.pro
;  Code determines the trajectory type for each severe stress event (MMM+2oC) that 
;   occurred during the analysis period (1985-2011) for each reef pixel on the GBR.
;   
; Developer: Scott F. Heron, scott.heron@noaa.gov
;======================================


pro AinsworthEtAl_2016_gbrReefs

  ; define temperature offset for future warming (=0 for historical analysis)
  TEMP_OFFSET = 0. ; [0., 0.25, 0.5, 1.0, 1.5, 2.0]

  path = '/'
  
  ; --------------------------------------------
  ;  CSV file "AinsworthEtAl_2016_SST_data.csv" contains LAT, LON, SST and MMM data.

  
  ; load GBR-wide SST time-series structures from proprietary data format (IDL).
  restore, path +'AinsworthEtAl_2016_PFV52_50km_GBRreefs.sav', restored_objects=vars
  i = ts_50km.pix_lon*2 +360
  j = ts_50km.pix_lat*2 +160  
  n_ts = n_elements(ts_50km.lon)
  
  ; read PFV52 MMM climatology
  pfr = PFV52_50km(/debug)
  mmm = pfr.clim(mm=mm)
  ; --------------------------------------------

  
  ; output text file
  openw, lun, path+'GBRreefs_above_red-offset'+$
    string(TEMP_OFFSET, format='(f0.2)') +'.txt', /get_lun
  printf, lun, "LAT,LON,YEAR,TYPE,b0,b1,b2,b3,b4,r5,r6,maxDHW,pulse_max,peak_max"

  n_noH = 0
  n_LH = 0
  n_LMLH = 0
  n_LHLH = 0

  for ct_ts = 0, n_ts-1 do begin

    lower = mmm[i[ct_ts], j[ct_ts]]
    bleaching_threshold = lower +1
    upper = lower +2
    
    loc_string = string(abs(ts_50km.lon[ct_ts]), format='(f05.1)')
    loc_string += (ts_50km.lon[ct_ts] lt 0) ? 'W_' : 'E_'
    loc_string += string(abs(ts_50km.lat[ct_ts]), format='(f04.1)')
    loc_string += (ts_50km.lat[ct_ts] lt 0) ? 'S' : 'N'
    print, loc_string, mmm[i[ct_ts], j[ct_ts]], lower, upper, $
      format='(a0, x, "MMM = ", f0.2, ", Lower = ", f0.2, ", Upper = ", f0.2)'

    ; consider each summer by the year in which summer ends (i.e., March)
    caldat, ts_50km.time, months, days, years
    start_year = 1985
    end_year = 2011; record ends in Dec2011, so last entire summer is 10/11; year[-1]
    prestressed = []
    prebigstressed = []
    unstressed = []
    bleaching_stressed = []
    stressed = []
    
    caldat, ts_50km.time, time_months, time_days, time_years
    
    for year = start_year, end_year do begin
      found = where( ts_50km.time ge julday(11, 1, year-1) and ts_50km.time le julday(5, 1, year), n_found )
      
      if n_found gt 0 then begin
        sst_cut = ts_50km.sst[ct_ts, found] +TEMP_OFFSET

        hs = sst_cut-lower
        hs[where(hs lt 1., /null)] = 0.
        dhw = hs
        for ct_dhw = 1, 23 do dhw[ct_dhw] += dhw[ct_dhw-1]
        for ct_dhw = 24, n_found-1 do dhw[ct_dhw] = total(hs[ct_dhw-23:ct_dhw])
        dhw *=0.5 ;data are twice-weekly, so scaling here required.
        
        found_lower = where( sst_cut ge lower, n_found_lower )
        found_upper = where( sst_cut ge upper, n_found_upper )
        
        ;Determine if prestress occurred and magnitude of prestress.
        prestress = -1
        prebigstress = -1
        if n_found_lower gt 1 then begin
          f = where( diff(found_lower) gt 1, nf)
          start_lower = found_lower[0]
          end_lower = found_lower[-1]
          if nf ge 1 then begin
            start_lower = [start_lower, found_lower[[1+where( diff(found_lower) gt 1 )]]]
            end_lower = [found_lower[[where( diff(found_lower) gt 1 )]], end_lower]
          endif
          n_lower = n_elements( start_lower )
          is_upper = bytarr( n_lower )
          
          for ct = 0, n_lower-1 do begin
            is_upper[ct] = total( found_upper ge start_lower[ct] and found_upper le end_lower[ct] ) gt 0
          endfor
          
          if n_lower gt 1 then begin
            if total(is_upper) ge 2 then prebigstress = found_upper[0] $
              else prestress = where( diff(is_upper) eq 1 )
            
          endif
        endif
        
        if prestress[0] ne -1 then prestressed = [prestressed, year]
        
        if found_upper[0] ne -1 then stressed = [stressed, year] else unstressed = [unstressed, year]
        
        if prebigstress[0] ne -1 then prebigstressed = [prebigstressed, year]
        
        found_bleaching = where( sst_cut gt bleaching_threshold )
        if found_bleaching[0] ne -1 then bleaching_stressed = [bleaching_stressed, year]
        
        ;typology of SST trajectory and extract times where lower and upper thresholds were crossed
        if n_found_upper gt 0 then begin
          if prestress eq -1 and prebigstress eq -1 then type='LH' else $
            if total(is_upper) eq 1 then type='LMLH' else type='LHLH'

            r5 = found_upper[0] ;, format='(C(CDI02, "-", CMoA, "-", CYI))')
            peak_max = max( sst_cut[r5:-1], r6i) -lower
            pulse_max = !values.f_nan
            r6 = found_upper[0]+r6i ;, format='(C(CDI02, "-", CMoA, "-", CYI))')
            
            if type eq 'LH' then begin
              b4 = string(ts_50km.time[found[found_lower[0]]], format='(C(CDI02, "-", CMoA, "-", CYI))')
              b3 = "-"
              b2 = "-"
              b1 = "-"
              b0 = "-"
            endif
            
            if type eq 'LMLH' then begin
              f = where( start_lower le found_upper[0] )
              b4 = start_lower[f[-1]]
              b0 = start_lower[f[-1]-1]
              pulse_max = max( sst_cut[b0:b4-1], index ) -lower
              b1 = b0 +index
              f = where( end_lower ge b1 )
              b2 = end_lower[f[0]]+1
              test = min( sst_cut[b2:b4], index )
              b3 = b2 +index
              
              ;if time between initial and maximum peaks (r6-b1) is >=68 days then re-type as LH.
              ; 68 days was determined by iteratively removing outliers of this interval, where outliers were
              ;  outside mean+2*stdv.  End result was mean=34.2, stdv=16.9 and mean+2*stdv=68.0.
              if ts_50km.time[found[r6]] -ts_50km.time[found[b1]] ge 68 then begin
                b3 = "-"
                b2 = "-"
                b1 = "-"
                b0 = "-"
                type = 'LH'
              endif else begin
                b0 = string(ts_50km.time[found[b0]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b1 = string(ts_50km.time[found[b1]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b2 = string(ts_50km.time[found[b2]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b3 = string(ts_50km.time[found[b3]], format='(C(CDI02, "-", CMoA, "-", CYI))')
              endelse
              b4 = string(ts_50km.time[found[b4]], format='(C(CDI02, "-", CMoA, "-", CYI))')
            endif
            
            if type eq 'LHLH' then begin
              r5 = found_upper[-1]
              index = -1
              while r5 ne found_upper[0] and r5-found_upper[index-1] eq 1 do begin
                r5 -=1
                index = where(found_upper eq r5)
              endwhile
              peak_max = max( sst_cut[r5:-1], r6i) -lower
              r6 = r5+r6i

              f = where( start_lower le r5 )
              b4 = start_lower[f[-1]]
              b0 = start_lower[f[-1]-1]
              pulse_max = max( sst_cut[b0:b4-1], index ) -lower
              b1 = b0 +index
              f = where( end_lower ge b1 )
              b2 = end_lower[f[0]]+1
              test = min( sst_cut[b2:b4], index )
              b3 = b2 +index

              ;if time between initial and maximum peaks (r6-b1) is >=60 days then re-type as LH.
              ; 60 days was determined by iteratively removing outliers of this interval, where outliers were
              ;  outside mean+2*stdv.  End result was mean=32.1, stdv=13.6 and mean+2*stdv=59.3.
              if ts_50km.time[found[r6]] -ts_50km.time[found[b1]] ge 60 then begin
                b3 = "-"
                b2 = "-"
                b1 = "-"
                b0 = "-"
                type = 'LH'
              endif else begin
                b0 = string(ts_50km.time[found[b0]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b1 = string(ts_50km.time[found[b1]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b2 = string(ts_50km.time[found[b2]], format='(C(CDI02, "-", CMoA, "-", CYI))')
                b3 = string(ts_50km.time[found[b3]], format='(C(CDI02, "-", CMoA, "-", CYI))')
              endelse
              b4 = string(ts_50km.time[found[b4]], format='(C(CDI02, "-", CMoA, "-", CYI))')
            endif
            r5 = string(ts_50km.time[found[r5]], format='(C(CDI02, "-", CMoA, "-", CYI))')
            r6 = string(ts_50km.time[found[r6]], format='(C(CDI02, "-", CMoA, "-", CYI))')

            printf, lun, ts_50km.lat[ct_ts], ts_50km.lon[ct_ts], year, type, $
              b0, b1, b2, b3, b4, r5, r6, max(dhw), pulse_max, peak_max, $
              format='(2(f0.1, ","), i0, 8(",", a0), ",", f0.1, 2(",", f0.2))'
            
            case 1 of
              type eq 'LH' : n_LH +=1
              type eq 'LMLH' : n_LMLH +=1
              type eq 'LHLH' : n_LHLH +=1
            endcase

        endif else begin
          n_noH +=1
          printf, lun, ts_50km.lat[ct_ts], ts_50km.lon[ct_ts], year, 'no_H', $
            '-', '-', '-', '-', '-', '-', '-', max(dhw), $
            format='(2(f0.1, ","), i0, 8(",", a0), ",", f0.1)'
        endelse
        
      endif
      
    endfor
    
    print, 'Stressed summers (at or above upper threshold = red line)'
    print, stressed
    print, 'Stressed summers that experienced a pre-stress above the blue line that relaxed and then exceeded the red line'
    print, prestressed
    print, ''
    
  endfor
  
  print, temp_offset, n_LH, n_LMLH, n_LHLH
  
  free_lun, lun
  
end; main