; Since NCL 6.2.0 ( April 2, 2014 ), the following scripts are automatically loaded
;***********************************************
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
;;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"
;***********************************************
begin
;************************************************
; read in data
;************************************************

  f1 = "Europe_bis_ncl.txt"

  nmdl   = 7                ; # of models
  EU     = readAsciiTable(f1,nmdl,"float",1)
  EU@_FillValue = -9999
  printVarSummary(EU)       ; [396] x [7]

;************************************************
; Create a 'year-month' (yyyymm) array for convenience;
; https://www.ncl.ucar.edu/Document/Functions/Contributed/yyyymm_time.shtml
;************************************************

  yrStrtTxt  = 1980
  yrLastTxt  = 2012
  yyyyTxt    = ispan(yrStrtTxt, yrLastTxt, 1) 
  print(yyyyTxt)               ; yyyyTxt(33)
    
  yyyymmTxt = yyyymm_time(yrStrtTxt, yrLastTxt, "integer")
  print(yyyymmTxt)             ; [time | 396]  ==> time: [198001..201212]

;************************************************
;  Desired subset yyyy and yyyymm
;************************************************

  yrStrt    = 2000
  yrLast    = 2012
  nyrs      = yrLast-yrStrt+1

  ymStrt    = yrStrt*100 +  1
  ymLast    = yrLast*100 + 12
  ntim      = nyrs*12

;************************************************
;  Convenience and to facilitate reordering via named dimensions
;************************************************

  dimEU     = dimsizes(EU)
  NTIM      = dimEU(0)         ; all times in file  (NTIM=396)
  NYRS      = NTIM/12          ; =dimsizes(yyyyTxt)

  EU!0      = "time"
  EU!1      = "model"
  EU&time   =  yyyymmTxt
  printVarSummary(EU)          ; [time | 396] x [model | 7]

;************************************************
;  "pretty print" the subset matrix; Use coordinate subscripting
;************************************************

  print(" ")
  printVarSummary(EU({ymStrt:ymLast}, :))
  print(" ")
  write_matrix (EU({ymStrt:ymLast}, :), nmdl+"f10.3", False)
  print(" ")

;************************************************
;  For each year, find the maximum NO2 for each model
;************************************************
  
  EUmax = new( (/NYRS,nmdl/), typeof(EU), EU@_FillValue)
  do nt=0,NTIM-1,12
    do nm=0,nmdl-1
       EUmax(nt/12,nm) = dim_max_n(EU(nt:nt+11,nm), 0)
    end do   ; nm
  end do     ; nt

  EUmax!0 = "yyyy" 
  EUmax!1 = "model"
  EUmax&yyyy = yyyyTxt
  print("")
  printVarSummary(EUmax)
  print("")
  print(yyyyTxt+"  "+EUmax(:,0)+"  "+EUmax(:,1)+"  "+EUmax(:,2) \
               +"  "+EUmax(:,3)+"  "+EUmax(:,4)+"  "+EUmax(:,5) \
               +"  "+EUmax(:,6)  )
  print("")

;************************************************
; plot multiple lines on same plot
; (a) place all into a single  mulidimensional array
;     or
; (b) use overlays
;************************************************
; Use (a) & (b):  (a) raw data, (b) overlay the mean 
;************************************************
; plotting parameters
;************************************************
  yrfrac = yyyymm_to_yyyyfrac(yyyymmTxt, 0.0)
  yrfrac!0    = "time"
  yrfrac&time = EU&time 
 ;printVarSummary(yrfrac)
  print(yrfrac)
  print("")

  wks    = gsn_open_wks ("png","Europe")                ; open workstation
  res              = True                   ; plot mods desired
  res@gsnDraw      = False                  ; don't draw
  res@gsnFrame     = False                  ; don't advance frame      
  res@tiMainString = "Multiple Models"      ; add title
  res@tiYAxisString= "Max NO2"              ; y-axis label
  res@vpHeightF    = 0.4                    ; change aspect ratio of plot
  res@vpWidthF     = 0.8                  
  res@vpXF         = 0.125                  ; start plot at x ndc coord
  res@tmXBFormat   = "f" 

 ;printVarSummary(yrfrac({ymStrt:ymLast}))  ; overview of x-axis being used
 ;printVarSummary(EU(model|:,{time|ymStrt:ymLast}))   ; overview of data being plotted
                                                      ; [model | 7] x [time | 156]
  plot   = gsn_csm_xy (wks,yrfrac({ymStrt:ymLast}) \
                          ,EU(model|:,{time|ymStrt:ymLast}),res) ; create plot

; compute average and add to plot

  EU_avg = dim_avg_n_Wrap(EU(model|:,{time|ymStrt:ymLast}),0)
  printVarSummary(EU_avg)

  plres = True
  plres@gsLineColor      = "blue"
  plres@gsLineThicknessF = 2.0
  stru = unique_string("EU")   ; will return a unique id
  plot@$stru$ = gsn_add_polyline(wks, plot, yrfrac({ymStrt:ymLast}),EU_avg, plres)
  
  draw(plot)
  frame(wks)   
end