#! /usr/bin/env nclsh
;
; $Id: plot_timeser 660 2022-09-22 10:06:49Z m221078 $
;

; ICON
;
; ---------------------------------------------------------------
; Copyright (C) 2004-2026, DWD, MPI-M, DKRZ, KIT, ETH, MeteoSwiss
; Contact information: icon-model.org
; See AUTHORS.TXT for a list of authors
; See LICENSES/ for license information
; SPDX-License-Identifier: BSD-3-Clause
; ---------------------------------------------------------------

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"
load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"

;
; Global settings
;

default_debug_mode = False
valid_modes = (/"monitoring", "presentation"/)
default_format = "pdf"
default_legend_mode = True
default_panel_rows = 4
default_output_file = "plot_timeser"
default_split_pages = False
default_units_pos = "top"
valid_units_pos = (/"top", "left", "none"/)
default_attribute_names = (/"code"/)
default_text_scale = 1.0
default_line_scale = 1.0
default_manifest = ""
default_time_format = "%Y-%N-%D~C~%H:%M:%S"
default_time_ticks = 33

debug_mode = default_debug_mode

hack_calendar_attribute_counter = 0

;
; Utility routines
;

;
; Write debug messages
;
undef("debug_print")
procedure debug_print(message)
begin
    if(debug_mode) then
        print("Debug: " + message)
    end if
end

;
; Get long name and units in a single string.
;
undef("get_long_name")
function get_long_name(data, default)
local result
begin
    result = default
    lname = get_valid_long_name_value(data)
    if(.not. ismissing(lname(0))) then
        result = lname(1)
    end if
    return result
end

undef("get_units")
function get_units(data, default)
local result, uname
begin
    result = default
    uname = get_valid_units_value(data)
    if(.not. ismissing(uname(0))) then
        result = uname(1)
    end if
    return result
end

;
; Check if a name belongs to a time coordinate
;
undef("is_time_name")
function is_time_name(datafile [1]: file, name)
begin
    return \
        isfilevardim(datafile, name, name) .and. \
        dimsizes(getfilevardimsizes(datafile, name)) .eq. 1 .and. \
        ( \
            ( isfilevaratt(datafile, name, "axis") .and. \
              datafile->$name$@axis .eq. "T" ) \
        .or. \
            ( isfilevaratt(datafile, name, "units") .and. \
              .not. ismissing( \
                  str_index_of_substr(datafile->$name$@units, " since ", 1) \
              ) ) \
        )
end

;
; Get name of time coordinate
;
undef("get_time_name")
function get_time_name(datafile [1]: file)
local names, name, var_max, var_idx
begin
  names = getfilevarnames(datafile)
  var_max = dimsizes(names)-1

  ; Coordinate vars have the same name for data and dimension, and are 1D.
  ; CF conforming time coords have an 'axis' attribute of 'T' or
  ; a relative time description in their 'units' attribute
  ; Just take the name of the first match.
  do var_idx = 0, var_max
    name = names(var_idx)
    if(is_time_name(datafile, name)) then
        return(name)
    end if
  end do
  name = default_fillvalue("string")
  name@_FillValue = name
  return(name)
end

;
; Create an empty list with default settings.
;
undef("empty_list")
function empty_list()
local result
begin
  result = NewList("fifo")
  ListSetType(result, "join") ; convert to an array easily (needed since 6.1)
  return(result)
end ; empty_list

;
; Remove in-sequence duplicates from list.
;
undef("uniq")
function uniq(xs[*])
local size, result, result_list, i
begin
    size = dimsizes(xs)
    result_list = empty_list()
    ListPush(result_list, ""+xs(size-1)) ; work-around for bug in ListPush
    do i = size-1, 0, 1
        if(xs(i) .ne. result_list[0]) then
            ListPush(result_list, ""+xs(i)) ; work-around for bug in ListPush
        end if
    end do
    result = result_list[:]
    return(result)
end ; uniq

;
; Simplify file names to non-constant parts for labeling.
;
undef("get_labels_from_file_names")
function get_labels_from_file_names(files[*]:string)
local \
    split_list, parts_num, file_max, file_idx, parts, part, part_index, \
    part_idx, first_pattern, mismatch, pattern, labels
begin

  split_list = "/\_.:"

  file_max = dimsizes(files)-1

  ; Singleton list is just passed back
  if(file_max .eq. 0) then
    files(file_max) = str_strip(str_join(str_split(files(file_max), split_list), " "))
    return(files)
  end if

  ; Split names at some of the usual separators, get size
  parts_num = 0
  do file_idx = 0, file_max
    parts_num = max((/parts_num, dimsizes(str_split(files(file_idx), split_list))/))
  end do
  ; ; Assume last part to be file extension
  ; parts_num = parts_num - 1

  ; Split names at some of the usual separators, store
  parts = new((/file_max+1, parts_num/), string)
  parts = ""
  do file_idx = 0, file_max
    part = str_split(files(file_idx), split_list)
    ; Put parts at end of array
    parts(file_idx,parts_num-dimsizes(part):) = part
    delete(part)
  end do

  ; Match up parts, taking all that are not found in all names
  part_index = empty_list()
  do part_idx = parts_num-1, 0, 1 ; reverse loop
    first_pattern = parts(0,part_idx)
    mismatch = False
    do file_idx = 1, file_max
        pattern = parts(file_idx,part_idx)
        ; if(first_pattern .eq. "" .and. pattern .ne. "") then
        ;     first_pattern = pattern
        ; end if
        ; if(pattern .ne. "" .and. first_pattern .ne. pattern) then
        if(first_pattern .ne. pattern) then
            mismatch = True
        end if
    end do
    if(mismatch) then
        ListPush(part_index, 0+part_idx) ; work-around for bug in ListPush
    end if
  end do

  ; Create labels from match list
  labels = new(file_max+1, string)
  do file_idx = 0, file_max
    labels(file_idx) = str_strip(str_join(uniq(parts(file_idx,part_index[:])), " "))
    if(labels(file_idx) .eq. "") then
        labels(file_idx) = "[baseline]"
    end if
  end do

  return(labels)

end ; get_labels_from_file_names

;
; Duplicate first element in list if list is a singleton
;
undef("fill2")
function fill2(xs[*])
local result
begin
    if(dimsizes(xs) .eq. 1) then
        result = (/ xs, xs /)
        copy_VarAtts(xs, result)
        result@filled = True
    else
        result = xs
        result@filled = False
    end if
    return(result)
end ; fill2

;
; Evaluate 'with' command line parameter
;
undef("push_with_parameter")
procedure push_with_parameter(with, file_list, label_list)
local size
begin
    size = dimsizes(with)
    if(size .gt. 1) then
        ListPush(file_list, ""+with(0))
        ListPush(label_list, str_join(with(1:size-1),","))
    else
        ListPush(file_list, ""+with)
        ListPush(label_list, ""+with)
    end if
end ; push_with_parameter

;
; Compare NCL versions
;
undef("version_cmp")
function version_cmp(version2)
local elem1, elem2, idx, temp
begin
  elem1 = stringtoint(str_split(get_ncl_version(), "."))
  elem2 = stringtoint(str_split(version2, "."))
  do idx = 0, (dimsizes(elem1) < dimsizes(elem2)) - 1
      temp = elem1(idx) - elem2(idx)
      if(temp .ne. 0) then
          return temp
      end if
  end do
  return 0
end ; version_cmp

undef("version_lt")
function version_lt(version2)
begin
  return version_cmp(version2) .lt. 0
end ; version_lt

undef("version_le")
function version_le(version2)
begin
  return version_cmp(version2) .le. 0
end ; version_le

undef("version_ge")
function version_ge(version2)
begin
  return version_cmp(version2) .ge. 0
end ; version_ge

;
; Work around for NCL bug in calendar routines
;
undef("hack_calendar_attribute")
procedure hack_calendar_attribute(time)
begin
  if(version_ge("6.1.0")) then
    if(isatt(time, "calendar") .and. time@calendar .eq. "proleptic_gregorian") then
        if(hack_calendar_attribute_counter .eq. 0) then
            print("Hey: replacing Proleptic Gregorian calendar by Gregorian calendar")
        end if
        if(hack_calendar_attribute_counter .eq. 1) then
            print("Hey: suppressing further Proleptic Gregorian calendar warnings")
        end if
        time@calendar = "gregorian"
        hack_calendar_attribute_counter = hack_calendar_attribute_counter + 1
    end if
  end if
end ; hack_calendar_attribute

;
; Main routine
;

begin

  ; Check command line arguments

  ; Check debug settings
  if(.not. isvar("debug")) then
    debug_mode = default_debug_mode
  else
    debug_mode = debug
  end if

  ; Evaluate mode
  if(isvar("mode")) then
    if(mode .eq. "presentation") then
        default_text_scale = 1.6
        default_line_scale = 7.0
        default_panel_rows = 1
        default_attribute_names = (/""/)
    else
    if(mode .eq. "monitoring") then
        default_text_scale = 1.1
        default_line_scale = 2.0
        default_panel_rows = 1
    else
        print("Oops: invalid mode '" + mode + "', must be one of " + \
              "'" + str_join(valid_modes, "', '") + "'")
        status_exit(1)
    end if
    end if
  end if

  ; Supply default for workstation format
  if(.not. isvar("format")) then
    format = default_format
  end if

  ; Supply default for panel rows
  if(.not. isvar("rows")) then
    panel_rows = default_panel_rows
  else
    panel_rows = rows
    delete(rows)
  end if

  ; Supply default for output file name
  if(.not. isvar("output")) then
    output_file = default_output_file
  else
    output_file = output
    delete(output)
  end if

  ; Supply default splitting of pages
  if(.not. isvar("split_pages")) then
    split_pages = default_split_pages
  end if

  ; Supply default for position of units' string
  if(.not. isvar("units")) then
    units_pos = default_units_pos
  else
    if(any(valid_units_pos .eq. units)) then
        units_pos = units
        delete(units)
    else
        print("Oops: invalid units position '" + units + "', must be one of " + \
              "'" + str_join(valid_units_pos, "', '") + "'")
        status_exit(1)
    end if
  end if

  ; Supply default for name of additional attribute
  if(.not. isvar("attr")) then
    attribute_names = default_attribute_names
  else
    attribute_names = attr
    delete(attr)
  end if

  ; Supply default text scale
  if(.not. isvar("text_scale")) then
    text_scale = default_text_scale
  end if

  ; Supply default line scale
  if(.not. isvar("line_scale")) then
    line_scale = default_line_scale
  end if

  ; Supply default manifest name
  if(.not. isvar("manifest")) then
    manifest = default_manifest
  end if

  ; Supply format for time labels
  if(.not. isvar("time_format")) then
    time_format = default_time_format
  end if
  if(.not. isvar("time_ticks")) then
    time_ticks = default_time_ticks
  end if

  ; Set distinct colors for consecutive line plots
  ;
  ; Specs legend:
  ;  0: Background (white)
  ;  1: Foreground (black)
  ; Saturated (a) colors:
  ;  2: Color  1a (red)
  ;  3: Color  2a (cyan)
  ;  4: Color  3a (yellow)
  ;  5: Color  4a (blue)
  ;  6: Color  5a (green)
  ;  7: Color  6a (magenta)
  ; 14: Color  7a (orange)
  ; 16: Color  8a (azure)
  ; 18: Color  9a (lime)
  ; 20: Color 10a (violet)
  ; 22: Color 11a (mint)
  ; 24: Color 12a (rose)
  ; Desaturated colors:
  ;  8: Color  1 (red)
  ;  9: Color  2 (cyan)
  ; 10: Color  3 (yellow)
  ; 11: Color  4 (blue)
  ; 12: Color  5 (green)
  ; 13: Color  6 (magenta)
  ; 15: Color  7 (orange)
  ; 17: Color  8 (azure)
  ; 19: Color  9 (lime)
  ; 21: Color 10 (violet)
  ; 23: Color 11 (mint)
  ; 25: Color 12 (rose)
  color_specs = (/ (/1.00, 1.00, 1.00/), \
                   (/0.00, 0.00, 0.00/), \
                                      \
                   (/1.00, 0.00, 0.00/), \
                   (/0.00, 0.80, 0.90/), \
                   (/1.00, 0.80, 0.00/), \
                   (/0.00, 0.00, 0.90/), \
                   (/0.00, 0.80, 0.00/), \
                   (/1.00, 0.00, 0.90/), \
                                         \
                   (/1.00, 0.50, 0.50/), \
                   (/0.50, 0.80, 0.90/), \
                   (/1.00, 0.80, 0.50/), \
                   (/0.50, 0.50, 0.90/), \
                   (/0.50, 0.80, 0.50/), \
                   (/1.00, 0.50, 0.80/), \
                                         \
                   (/1.00, 0.50, 0.00/), \
                   (/1.00, 0.65, 0.50/), \
                                         \
                   (/0.00, 0.50, 0.90/), \
                   (/0.50, 0.65, 0.90/), \
                                         \
                   (/0.50, 0.80, 0.00/), \
                   (/0.75, 0.80, 0.50/), \
                                         \
                   (/0.50, 0.00, 0.90/), \
                   (/0.75, 0.50, 0.90/), \
                                         \
                   (/0.00, 0.80, 0.50/), \
                   (/0.50, 0.80, 0.70/), \
                                         \
                   (/1.00, 0.00, 0.50/), \
                   (/1.00, 0.50, 0.70/), \
                                         \
                   (/0,0,0/) /) ; end-of-list dummy

  ; Explicitly assign files to colors
  file_list = empty_list()
  label_list = empty_list()
  color_list = empty_list()
  if(isvar("with1a")) then
    push_with_parameter(with1a, file_list, label_list)
    ListPush(color_list, 2)
  end if
  if(isvar("with2a")) then
    push_with_parameter(with2a, file_list, label_list)
    ListPush(color_list, 3)
  end if
  if(isvar("with3a")) then
    push_with_parameter(with3a, file_list, label_list)
    ListPush(color_list, 4)
  end if
  if(isvar("with4a")) then
    push_with_parameter(with4a, file_list, label_list)
    ListPush(color_list, 5)
  end if
  if(isvar("with5a")) then
    push_with_parameter(with5a, file_list, label_list)
    ListPush(color_list, 6)
  end if
  if(isvar("with6a")) then
    push_with_parameter(with6a, file_list, label_list)
    ListPush(color_list, 7)
  end if
  if(isvar("with7a")) then
    push_with_parameter(with7a, file_list, label_list)
    ListPush(color_list, 14)
  end if
  if(isvar("with8a")) then
    push_with_parameter(with8a, file_list, label_list)
    ListPush(color_list, 16)
  end if
  if(isvar("with9a")) then
    push_with_parameter(with9a, file_list, label_list)
    ListPush(color_list, 18)
  end if
  if(isvar("with10a")) then
    push_with_parameter(with10a, file_list, label_list)
    ListPush(color_list, 20)
  end if
  if(isvar("with11a")) then
    push_with_parameter(with11a, file_list, label_list)
    ListPush(color_list, 22)
  end if
  if(isvar("with12a")) then
    push_with_parameter(with12a, file_list, label_list)
    ListPush(color_list, 24)
  end if
  if(isvar("with1")) then
    push_with_parameter(with1, file_list, label_list)
    ListPush(color_list, 8)
  end if
  if(isvar("with2")) then
    push_with_parameter(with2, file_list, label_list)
    ListPush(color_list, 9)
  end if
  if(isvar("with3")) then
    push_with_parameter(with3, file_list, label_list)
    ListPush(color_list, 10)
  end if
  if(isvar("with4")) then
    push_with_parameter(with4, file_list, label_list)
    ListPush(color_list, 11)
  end if
  if(isvar("with5")) then
    push_with_parameter(with5, file_list, label_list)
    ListPush(color_list, 12)
  end if
  if(isvar("with6")) then
    push_with_parameter(with6, file_list, label_list)
    ListPush(color_list, 13)
  end if
  if(isvar("with7")) then
    push_with_parameter(with7, file_list, label_list)
    ListPush(color_list, 15)
  end if
  if(isvar("with8")) then
    push_with_parameter(with8, file_list, label_list)
    ListPush(color_list, 17)
  end if
  if(isvar("with9")) then
    push_with_parameter(with9, file_list, label_list)
    ListPush(color_list, 19)
  end if
  if(isvar("with10")) then
    push_with_parameter(with10, file_list, label_list)
    ListPush(color_list, 21)
  end if
  if(isvar("with11")) then
    push_with_parameter(with11, file_list, label_list)
    ListPush(color_list, 23)
  end if
  if(isvar("with12")) then
    push_with_parameter(with12, file_list, label_list)
    ListPush(color_list, 25)
  end if
  if(isvar("with0")) then
    push_with_parameter(with0, file_list, label_list)
    ListPush(color_list, 1)
  end if

  ; Filenames are mandatory
  if(ListCount(file_list) .eq. 0) then
    print("Invalid number of parameters")
    status_exit(1)
  end if
  files = file_list[:]
  labels = label_list[:]
  colors = color_list[:]

  ; Supply default for legend mode, depending of number of files
  if(.not. isvar("legend")) then
    legend_mode = dimsizes(files) .gt. 1 .and. default_legend_mode
  else
    legend_mode = legend
    delete(legend)
  end if

  ; Open input files

  ; Read files
  datafiles = addfiles(files, "r")
  file_max = dimsizes(files)-1

  if(legend_mode) then
    ; Simplify file names to non-constant parts for labeling
    labels = get_labels_from_file_names(labels)
  end if

  ; Get time coordinate name from first file,
  ; assuming same variable structure across files.
  ; Read in all time axis values and get min and max

  time_name = get_time_name(datafiles[file_max])
  if(ismissing(time_name)) then
    print("Sorry: cannot find valid time coordinate in '" + files(file_max) + "'")
    status_exit(1)
  end if

  ; Traverse files back to front as list only provides 'push' operation

  debug_print("scan time axes")
  debug_print("file '" + files(file_max) + "'")
  time = datafiles[file_max]->$time_name$
  hack_calendar_attribute(time)
  time_units = time@units
  x = [/ fill2(time) /]
  x_min = min(x[0])
  x_max = max(x[0])
  debug_print("range [" + x_min + ", " + x_max + "]")
  delete(time)
  if(file_max .gt. 0) then
    do file_idx = file_max-1, 0, 1
        debug_print("file '" + files(file_idx) + "'")
        time = datafiles[file_idx]->$time_name$
        hack_calendar_attribute(time)
        ListPush(x, fill2(cd_convert(time, time_units)))
        x_min = min((/x_min, min(x[0])/))
        x_max = max((/x_max, max(x[0])/))
        debug_print("range [" + x_min + ", " + x_max + "]")
        delete(time)
    end do
  end if
  x_values := fspan(x_min, x_max, time_ticks)
  x_values@units = time_units

  ; Open graphic workstation

  if(panel_rows .eq. 1) then
    format@wkPaperHeightF = 8.3
    format@wkPaperWidthF = 11.7
  else
    format@wkPaperSize = "A4"
  end if
  if( version_lt("6.1.0") .and. any((/"ps", "pdf"/) .eq. format) .or. \
      any((/"oldps", "oldpdf"/) .eq. format) ) \
  then
      format@wkFullBackground = True
  end if

  ; Setup resources for every plot

  res = True

  res@gsnFrame = False
  res@gsnDraw = False

  if(panel_rows .eq. 1) then
    res@vpHeightF = 0.6
    res@vpWidthF = 0.8
  else
    res@vpHeightF = 0.8
    res@vpWidthF = 0.6
  end if
  res@vpHeightF = res@vpHeightF / panel_rows
  res@gsnStringFontHeightF = 0.012 * text_scale

  res@tmBorderThicknessF = 0.5 * line_scale
  res@tmXBMajorThicknessF = res@tmBorderThicknessF
  res@tmYLMajorThicknessF = res@tmBorderThicknessF
  res@tmXBMinorThicknessF = 0.75 * res@tmBorderThicknessF
  res@tmYLMinorThicknessF = 0.75 * res@tmBorderThicknessF

  res@tmXMajorGrid = True
  res@tmXMajorGridLineDashPattern = 2
  res@tmXMajorGridThicknessF = res@tmBorderThicknessF

  res@tmYMajorGrid = True
  res@tmYMajorGridLineDashPattern = res@tmXMajorGridLineDashPattern
  res@tmYMajorGridThicknessF = res@tmBorderThicknessF

  res@tmEqualizeXYSizes = True
  res@tmXBLabelFontHeightF = res@gsnStringFontHeightF

  res@tmLabelAutoStride = True

  res@trXMinF = tofloat(x_min)
  res@trXMaxF = tofloat(x_max)

  ; Setup time labels

  res_time = True
  res_time@ttmFormat = time_format
  time_axis_labels(x_values, res, res_time)

  ; Setup page layout (panel)

  panel_res = True
  panel_res@gsnMaximize      = True
  panel_res@gsnPaperOrientation = "portrait"
  ;;; panel_res@gsnPanelYWhiteSpacePercent = 5

  plot_row = 0
  plots = new(panel_rows, graphic)

  if(legend_mode) then

    ; Setup panel legend (actually a list of text items)

    tx_res = True

    ; Set font size
    tx_res@txFontHeightF = 0.875 * res@gsnStringFontHeightF

    ; Allocate text items
    texts = new(file_max+1, graphic)

  end if

  ; Read variable names from first file,
  ; assuming same variable structure across files.

  names = getfilevarnames(datafiles[file_max])
  var_max = dimsizes(names)-1

  ; Initialize manifest for pages
  manifest_list = empty_list()

  do var_idx = 0, var_max

    debug_print("scan '" + names(var_idx) + "' axes")
    var_dims = getfilevardimsizes(datafiles[file_max], names(var_idx))
    debug_print("dims(" + names(var_idx) + ") = " + var_dims)
    var_rank = dimsizes(var_dims)
    debug_print("is not string data: " + \
        (getfilevartypes(datafiles[file_max], names(var_idx)) .ne. "string") )
    debug_print("has dim '" + time_name + "': " + \
        (isdim(datafiles[file_max]->$names(var_idx)$, time_name)) )
    debug_print("has not time name: " + \
        (.not. is_time_name(datafiles[file_max], names(var_idx))) )
    debug_print("has rank 1: " + \
        (product(var_dims) .eq. max(var_dims)) )
    if( getfilevartypes(datafiles[file_max], names(var_idx)) .ne. "string" .and.\
        isdim(datafiles[file_max]->$names(var_idx)$, time_name) .and. \
        .not. is_time_name(datafiles[file_max], names(var_idx)) .and. \
        product(var_dims) .eq. max(var_dims) ) \ ; i.e. effectively rank 1
    then

        ; Read in all variable values and get min and max

        ; Read files back to front as list only provides 'push' operation
        y_type = getfilevartypes(datafiles[file_max], names(var_idx))
        y_fill = getFillValue(datafiles[file_max]->$names(var_idx)$)
        y_min = new(1, y_type, y_fill)
        y_max = y_min
        y_avg = y_max
        y = empty_list()
        do file_idx = file_max, file_max, 1
            if(isfilevar(datafiles[file_idx], names(var_idx))) then
                debug_print("file '" + files(file_idx) + "'")
                y_tmp = ndtooned(datafiles[file_idx]->$names(var_idx)$)
                copy_VarAtts(datafiles[file_idx]->$names(var_idx)$, y_tmp)
                ListPush(y, fill2(y_tmp))
                delete(y_tmp)
                y_min = min(y[0])
                y_max = max(y[0])
                y_avg = avg(y[0])
                debug_print("range [" + y_min + ", " + y_max + "]")
            else
                ; Skip missing variables
                ListPush(y, fill2(new(dimsizes(x[file_idx]), y_type, y_fill)))
            end if
        end do
        do file_idx = file_max, 0, 1
            if(isfilevar(datafiles[file_idx], names(var_idx))) then
                debug_print("file '" + files(file_idx) + "'")
                y_tmp = ndtooned(datafiles[file_idx]->$names(var_idx)$)
                copy_VarAtts(datafiles[file_idx]->$names(var_idx)$, y_tmp)
                ListPush(y, fill2(y_tmp))
                delete(y_tmp)
                y_min = min((/y_min, min(y[0])/))
                y_max = max((/y_max, max(y[0])/))
                y_avg = avg((/y_avg, avg(y[0])/))
                debug_print("range [" + y_min + ", " + y_max + "]")
            else
                ; Skip missing variables
                ListPush(y, fill2(new(dimsizes(x[file_idx]), y_type, y_fill)))
            end if
        end do
        delete(y_fill)
        delete(y_type)

        ; Beautify x scaling
        if(x_max .eq. x_min) then
            res@trXMinF = tofloat(x_min - 0.1)
            res@trXMaxF = tofloat(x_max + 0.1)
        end if

        ; Beautify y scaling
        if(y_max .eq. y_min) then
            debug_print("y range is empty")
            if(y_max .eq. 0.) then
                res@trYMinF = -0.01
                res@trYMaxF = 0.01
            else
                res@trYMinF = tofloat(y_min - 0.01*abs(y_min))
                res@trYMaxF = tofloat(y_max + 0.01*abs(y_max))
            end if
        else
            debug_print("y range spans " + (y_max-y_min))
            if((y_max-y_min) < 0.0001*abs(y_min) .or. \
               (y_max-y_min) < 0.0001*abs(y_max)) \
            then
                res@trYMinF = tofloat(y_min - 0.0001*abs(y_min))
                res@trYMaxF = tofloat(y_max + 0.0001*abs(y_max))
            else
                res@trYMinF = tofloat(y_min - 0.01*(y_max-y_min))
                res@trYMaxF = tofloat(y_max + 0.01*(y_max-y_min))
            end if
        end if
        debug_print("y range adjusted to [" + res@trYMinF + ", " + res@trYMaxF + "]")

        long_name = get_long_name(y[file_max], names(var_idx))
        units = get_units(y[file_max], "")

        if(units_pos .eq. "top") then
            if(str_is_blank(units)) then
                res@gsnLeftString = long_name
            else
                res@gsnLeftString = long_name + " [" + units + "]"
            end if
            res@tiYAxisString = ""
        else
            if(units_pos .eq. "left") then
                res@gsnLeftString = long_name
                res@tiYAxisString = units
                res@tiYAxisFontHeightF = res@gsnStringFontHeightF
            else
                res@gsnLeftString = long_name
                res@tiYAxisString = ""
            end if
        end if

        attribute = ""
        res@gsnRightString = attribute;
        att_max = dimsizes(attribute_names) - 1
        do att_idx = 0, att_max
            attribute_name = attribute_names(att_idx)
            if(isatt(y[file_max], attribute_name)) then
                if(str_is_blank(attribute)) then
                    attribute = "" + y[file_max]@$attribute_name$
                else
                    attribute = attribute + ", " + y[file_max]@$attribute_name$
                end if
            end if
        end do
        if(.not. str_is_blank(attribute)) then
            res@gsnRightString = "(" + attribute + ")"
        end if

        file_idx = 0
        num_defined = num(.not. ismissing(y[file_idx]))
        debug_print("file "+ file_idx + " has " + num_defined + " defined values")
        if(num_defined .eq. 2) then
            res@xyMarkLineMode = "Markers"
            res@xyMarkerColor = colors(file_idx)
        else
            res@xyMarkLineMode = "Lines"
            res@xyLineColor = colors(file_idx)
        end if

        if(colors(file_idx) .eq. 1) then
            res@xyLineThicknessF = 0.5 * line_scale
        else
            res@xyLineThicknessF = 1 * line_scale
        end if

        ; @todo Resolve overlay issue with reference lines
        ; if(.not. ismissing(y_avg)) then
        ;     res@gsnYRefLineColor = 31
        ;     res@gsnYRefLine = y_avg
        ; else
        ;     if(isatt(res, "gsnYRefLineColor")) then
        ;         delete(res@gsnYRefLineColor)
        ;     end if
        ;     if(isatt(res, "gsnYRefLine")) then
        ;         delete(res@gsnYRefLine)
        ;     end if
        ; end if
        ; res@gsnYRefLine = avg(y[file_max])
        ; res@gsnYRefLineColor = res@xyLineColor

        ; Set up new page
        if(plot_row .eq. 0) then

            ; Append variable name for split pages
            if(split_pages) then
                base_name = output_file + "_" + names(var_idx)
            else
                base_name = output_file
            end if

            ; Open and initialize canvas
            if(.not. isvar("wks")) then

                wks = gsn_open_wks(format, base_name)
                gsn_define_colormap(wks, color_specs)
                if(legend_mode) then
                  ; Initialize text items using plot colors
                  do file_idx = 0, file_max
                      tx_res@txFontColor = colors(file_idx)
                      texts(file_idx) = gsn_create_text(wks, labels(file_idx), tx_res)
                  end do
                end if

            end if; wks exists

            ; Add info to manifest
            idx = ListCount(manifest_list)
            ListPush(manifest_list, base_name+":"+idx+":"+format+":"+ \
                     names(var_idx)+":"+attribute+":"+res@gsnLeftString)

        end if; first plot on page

        debug_print("plotting row " + plot_row + ", file 0")
        plots(plot_row) = gsn_csm_xy(wks, x[0], y[0], res)

        do file_idx = 1, file_max
            num_defined = num(.not. ismissing(y[file_idx]))
            debug_print("file "+ file_idx + " has " + num_defined + " defined values")
            if(num_defined .gt. 0) then

                if(num_defined .eq. 2) then
                    res@xyMarkLineMode = "Markers"
                    res@xyMarkerColor = colors(file_idx)
                else
                    res@xyMarkLineMode = "Lines"
                    res@xyLineColor = colors(file_idx)
                end if

                if(colors(file_idx) .eq. 1) then
                    res@xyLineThicknessF = 0.5 * line_scale
                else
                    res@xyLineThicknessF = 1 * line_scale
                end if

                ; @todo Resolve overlay issue with reference lines
                ; res@gsnYRefLine = avg(y[file_idx])
                ; res@gsnYRefLineColor = res@xyLineColor

                debug_print("overlaying row " + plot_row + ", file " + file_idx)
                overlay(plots(plot_row), \
                        gsn_csm_xy(wks, x[file_idx], y[file_idx], res))

            end if
        end do

        delete(y)
        delete(y_avg)
        delete(y_max)
        delete(y_min)

        plot_row = plot_row + 1

    end if ; is valid variable with effective rank 1
    delete(var_dims) ; Delete because next variable might be different

    ; Do actual output at end of page or file

    if((plot_row .eq. panel_rows .or. var_idx .eq. var_max) .and. \
       .not. ismissing(plots(0)) .and. isvar("wks")) \
    then

        if(legend_mode) then
            ; Add legend texts as annotation into upper left corner of first plot
            am_res = True
            am_res@amJust = "TopLeft"
            am_res@amParallelPosF = -0.5 + tx_res@txFontHeightF
            annos = new(file_max+1, graphic)
            do file_idx = 0, file_max
                am_res@amOrthogonalPosF = -0.5 + (tx_res@txFontHeightF + \
                    file_idx*1.5/res@vpHeightF/panel_rows*tx_res@txFontHeightF)\
                    *panel_rows
                    ;;; file_idx*2.5*tx_res@txFontHeightF)*panel_rows
                annos(file_idx) = gsn_add_annotation(plots(0), texts(file_idx), am_res)
            end do; file_idx
        end if

        gsn_panel(wks, plots, (/panel_rows, 1/), panel_res)
        plots = getFillValue(plots)
        plot_row = 0

        ; Close current workstation to end page if splitting is requested
        if(split_pages) then
            delete(wks)
        end if

    end if

  end do; var_idx

  ; Write manifest for pages
  if(manifest .ne. "") then
    count = ListCount(manifest_list)
    asciiwrite(manifest, manifest_list[count-1:0])
  end if


end
