Changed ut_ to cd_ - Wed May 11 00:00:14 EDT 2016

econ/nca_spot.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; nca_spot.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; May 2016
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Draw economic plots for BAMS article
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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/csm/shea_util.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+load "sub.draw_spot.ncl"
+load "sub.draw_storage.ncl"
+load "sub.draw_nca.ncl"
+
+
+begin ; main
+  print_clock( "Here we go!" )
+
+  ; These are some parameters that could be useful to have up top
+
+  plotType = "png"
+  plotDpi  = 200
+
+  plotName = "figures/nca_spot"
+
+  fontHeightF = 0.02
+
+  timeUnits = "days since 1800-01-01 00:00:00"
+  minTime = cd_inv_calendar( 2011, 04, 01, 00, 0, 0, timeUnits, 0 )
+  maxTime = cd_inv_calendar( 2016, 04, 01, 00, 0, 0, timeUnits, 0 )
+
+  mwePath = "/home/carl/data/ghcn-daily/jared/results_NCA-US.nc"
+  mweVar  = "tavgAnom"
+
+  ; Customize panel
+  panRes = True
+  panRes@gsnPanelBottom = 0.10
+; panRes@gsnPanelYWhiteSpacePercent = 5
+
+  print_clock( "Drawing the plot" )
+
+  ; ...allows png or gif to work
+  if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
+    plotTypeLocal = "eps"
+  else
+    plotTypeLocal = plotType 
+  end if
+
+  ; ...open the workstation
+  wks = gsn_open_wks( plotTypeLocal, plotName )
+
+  print_clock( "Drawing the MWE plot" )
+  mwePlot = draw_nca( wks, minTime, maxTime, mwePath, mweVar, 2 )
+
+  print_clock( "Drawing the price plot" )
+  pricePlot = draw_spot( wks, minTime, maxTime, fontHeightF )
+; gsn_panel( wks, (/ mwePlot, pricePlot /), (/ 2, 1 /), panRes )
+
+  draw(mwePlot)
+  draw(pricePlot)
+  frame(wks)
+
+  print_clock( "Convert the image, if necessary" )
+  delete(wks)
+  if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
+    system( "convert -trim -border 5x5 -bordercolor white " \
+        + "+repage -density " + plotDpi + " " \\
+        + plotName + ".eps " + plotName + "." + plotType )
+  end if
+
+
+end; main
+

econ/storage.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; storage.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; March 2014
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Draw economic plots for BAMS article
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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/csm/shea_util.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+load "sub.draw_storage.ncl"
+
+
+begin ; main
+  print_clock( "Here we go!" )
+
+  ; These are some parameters that could be useful to have up top
+
+  plotType = "png"
+  plotDpi  = 200
+
+  plotName = "figures/storage"
+
+  fontHeightF = 0.02
+  allYears = ispan( 1994, 2015, 1 )
+  specialYears = (/ 2012, 2013, 2014, 2015, 2016 /)
+  startDDD = 244
+  startDDD = 152
+
+
+  print_clock( "Drawing the plot" )
+
+  ; ...allows png or gif to work
+  if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
+    plotTypeLocal = "eps"
+  else
+    plotTypeLocal = plotType 
+  end if
+
+  ; ...open the workstation
+  wks = gsn_open_wks( plotTypeLocal, plotName )
+
+  print_clock( "Drawing the background storage plot" )
+  allStoragePlot = draw_storage( wks, fontHeightF, allYears, startDDD )
+  print_clock( "Drawing the special storage plot" )
+  specialStoragePlot = draw_storage( wks, fontHeightF, specialYears, \
+      startDDD )
+; overlay( allStoragePlot, specialStoragePlot )
+  draw(allStoragePlot)
+  frame(wks)
+
+
+  print_clock( "Convert the image, if necessary" )
+  delete(wks)
+  if( isStrSubset( plotType, "png" ).or.isStrSubset( plotType, "gif" ) ) then
+    system( "convert -trim -border 5x5 -bordercolor white " \
+        + "+repage -density " + plotDpi + " " \\
+        + plotName + ".eps " + plotName + "." + plotType )
+    if( .not.isStrSubset( plotType, "e" ) ) then
+      system( "rm -f " + plotName + ".eps" )
+    end if
+  end if
+
+
+end; main
+

econ/sub.draw_nca.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; sub.draw_nca.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; March 2014
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Draw bar chart of a nca time series
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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/csm/shea_util.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+load "$CJS_NCL_LIB/lib.cjs_graphics.ncl"
+
+
+undef ( "draw_nca" )
+function draw_nca( \
+  io_wks         [1] : graphic, \
+  i_minTime      [1] : numeric, \
+  i_maxTime      [1] : numeric, \
+  i_inPath       [1] : string , \
+  i_varName      [1] : string , \
+  i_nRegion      [1] : integer \
+)
+begin ; draw_nca
+  ; These are some parameters that could be useful to have up top
+  fontHeightF = 0.025
+
+  inFile = addfile( i_inPath, "r" )
+  data = inFile->$i_varName$(i_nRegion,{i_minTime:i_maxTime})
+  data = runave( data, 30, 0 )
+
+  ; ...set up labels on the time axis
+  tickMonths = ispan( 0, 500, 1 )
+  tickZeros  = tickMonths * 0
+  tickDays   = tickZeros + 1
+  tickYears  = toint( cd_string( i_minTime, "%Y" ) ) + tickMonths / 12
+  tickMonths = 1 + ( tickMonths % 12 )
+  tmXBMinorValues := cd_inv_calendar( tickYears, tickMonths, tickDays, \
+      tickZeros, tickZeros, tickZeros, i_minTime@units, 0 )
+
+  res = True
+  res@gsnLeftString = "Northeast Temperature and Natural Gas Prices"
+  res@tmYROn = False
+  res@trYMaxF = 8
+  res@gsnYRefLine = 0
+  res@xyLineColor = -1
+  res@tmXBLabelAngleF = 45
+  res@tmXBLabelJust = "CenterRight"
+  res@tiXAxisString = ""
+  res@vpHeightF = 0.4
+  res@trYMinF = -12
+  res@trYMaxF =  12
+  
+
+  resTick = True
+  resTick@ttmFormat      = "%c-%Y"
+  resTick@ttmAxis        = "XB"
+  ttmValues = new( (/ 50, 6 /), integer )
+  ttmValues(:,0)   = toint( cd_string( i_minTime, "%Y" ) ) \
+      + ispan( 0, 49, 1 ) / 2
+  ttmValues(:,1)   = toint( cd_string( i_minTime, "%N" ) ) \
+      + ispan( 0, 49, 1 ) * 6 % 12
+  ttmValues(:,1)   = where( ttmValues(:,1).gt.12, \
+      ttmValues(:,1) - 12, ttmValues(:,1) )
+  ttmValues(:,2)   = 1
+  ttmValues(:,3:)  = 0
+  resTick@ttmValues = ttmValues
+; resTick@ttmMajorStride = 366
+; resTick@ttmMinorStride = 122
+  time_axis_labels( data&time,  res, resTick )
+  res@tmXBMinorValues := tmXBMinorValues
+
+  print_clock( "Drawing the plot" )
+
+  printMinMax(data,True)
+  retVal = cjs_draw_xy( io_wks, data&time, data, res )
+
+  return(retVal)
+end; draw_nca
+

econ/sub.draw_spot.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; sub.draw_spot.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; February 2014
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Draw a time series of the Henry Hub spot prices
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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/csm/shea_util.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+
+
+undef ( "draw_spot" )
+function draw_spot( \
+  io_wks[1] : graphic, \
+  i_minTime[1] : numeric , \
+  i_maxTime[1] : numeric , \
+  i_fontHeightF[1] : numeric \
+)
+begin ; draw_spot
+
+  ; These are some parameters that could be useful to have up top
+; inPath = "/home/carl/data/econ/henry_hub_spot.nc"
+  inPath = "/home/carl/data/econ/nymex_futures.nc"
+  inFile = addfile( inPath, "r" )
+  data = inFile->price({i_minTime:i_maxTime})
+  data = runave( data, 30, 0 )
+
+
+  ; Customize base plot
+  res = True
+  res@gsnDraw  = False
+  res@gsnFrame = False
+; res@gsnLeftString  = "b) " + data@long_name
+; res@gsnLeftString  = "b) Natural Gas Futures Price"
+  res@gsnLeftString = ""
+  res@tiYAxisString = data@units
+  res@trXMinF = i_minTime
+  res@trXMaxF = i_maxTime
+  res@vpHeightF = 0.4
+
+  ; ...standardize the font sizes
+  res@gsnLeftStringFontHeightF  = i_fontHeightF
+  res@gsnRightStringFontHeightF = i_fontHeightF
+  res@tiMainFontHeightF         = i_fontHeightF * 1.3
+  res@tiYAxisFontHeightF        = i_fontHeightF
+  res@tiXAxisFontHeightF        = i_fontHeightF
+  res@tmXBLabelFontHeightF      = i_fontHeightF
+  res@tmYLLabelFontHeightF      = i_fontHeightF
+  res@tmYLOn = False
+  
+  res@gsnLeftString = "Northeast Temperature and Natural Gas Prices"
+
+  ; ...set up an x-y line graph
+  res@xyMarkLineMode      = "Lines"
+  res@xyMonoDashPattern   = True
+  res@xyLineColors        = "darkolivegreen4"
+  res@xyMarkerColors      = res@xyLineColors
+  res@xyMarkers           = (/  16,   6,  2,  5,  7,  8,  9, 12 /)
+  res@xyMarkerThicknesses = (/   1,   1,  1,  1,  1,  1,  1,  1 /) * 2
+  res@xyLineThicknesses   = (/   2,   2,  1,  1,  1,  1,  1,  1 /) * 2
+  res@xyMarkerSizes       = (/ 1.5, 1.5,  1,  1,  1,  1,  1,  1 /) * 0.007
+
+  res@tmXBLabelAngleF = 45
+  res@tmXBLabelJust = "CenterRight"
+
+  ; ...make gridlines look nice
+  res@tmYMajorGridLineDashPattern = 2.
+  res@tmYMajorGridThicknessF      = 1.
+
+  res@trYMinF = 0
+  res@trYMaxF = 6
+  res@tmYRLabelsOn = True
+  res@tiYAxisSide = "Right"
+
+  ; ...set up labels on the time axis
+  tickMonths = ispan( 0, 500, 1 )
+  tickZeros  = tickMonths * 0
+  tickDays   = tickZeros + 1
+  tickYears  = toint( cd_string( i_minTime, "%Y" ) ) + tickMonths / 12
+  tickMonths = 1 + ( tickMonths % 12 )
+  tmXBMinorValues := cd_inv_calendar( tickYears, tickMonths, tickDays, \
+      tickZeros, tickZeros, tickZeros, i_minTime@units, 0 )
+
+  resTick = True
+  resTick@ttmFormat      = "%c-%Y"
+  resTick@ttmAxis        = "XB"
+  ttmValues = new( (/ 50, 6 /), integer )
+  ttmValues(:,0)   = toint( cd_string( i_minTime, "%Y" ) ) \
+      + ispan( 0, 49, 1 ) / 2
+  ttmValues(:,1)   = toint( cd_string( i_minTime, "%N" ) ) \
+      + ispan( 0, 49, 1 ) * 6 % 12
+  ttmValues(:,1)   = where( ttmValues(:,1).gt.12, \
+      ttmValues(:,1) - 12, ttmValues(:,1) )
+  ttmValues(:,2)   = 1
+  ttmValues(:,3:)  = 0
+  resTick@ttmValues = ttmValues
+; resTick@ttmMajorStride = 366
+; resTick@ttmMinorStride = 122
+  time_axis_labels( data&time,  res, resTick )
+  res@tmXBMinorValues := tmXBMinorValues
+
+  printMinMax(data,True)
+  retVal = gsn_csm_xy( io_wks, data&time, data, res )
+
+  return(retVal)
+end; draw_spot
+

econ/sub.draw_storage.ncl

+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; sub.draw_storage.ncl
+; Carl Schreck (cjschrec@ncsu.edu)
+; March 2014
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Description: Empty file to start NCL scripts
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+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/csm/shea_util.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/wrf/WRFUserARW.ncl"
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/diagnostics_cam.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/time_axis_labels.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/contrib/cd_string.ncl"
+load "$CJS_NCL_LIB/print_clock.ncl"
+
+
+undef ( "draw_storage" )
+function draw_storage( \
+  io_wks        [1] : graphic, \
+  i_fontHeightF [1] : numeric, \
+  i_plotYears   [*] : integer, \
+  i_startDDD    [1] : integer \
+)
+begin ; draw_storage
+
+  maxNddd = 60
+  i_startDDD@units = "days since 1-1-1"
+  nPlotYears = dimsizes(i_plotYears)
+
+  ; Read the data
+  inPath = "/home/carl/data/econ/ng_storage.nc"
+  inFile = addfile( inPath, "r" )
+  data = inFile->storage
+  dataYear = toint( cd_string( data&time, "%Y" ) )
+  dataDDD  = toint( cd_string( data&time, "%J" ) )
+  dataYear = where( dataDDD.ge.i_startDDD, dataYear+1, dataYear )
+  dataDDD = where( dataDDD.ge.i_startDDD, dataDDD, dataDDD+365 )
+  allYears = ispan( min(dataYear), max(dataYear), 1 )
+
+  ; reorder the data
+  do y = 0, dimsizes(allYears)-1
+    originalInd := ind( dataYear.eq.allYears(y) ) 
+    midInd = maxind(dataDDD(originalInd))
+    if( midInd.ne.dimsizes(originalInd)-1 ) then
+      springInd := originalInd(midInd+1:)
+      fallInd := originalInd(:midInd)
+      reorderedInd := originalInd
+      reorderedInd(:dimsizes(springInd)-1) = springInd
+      reorderedInd(dimsizes(springInd):) = fallInd
+      data(originalInd) = data(reorderedInd)
+      dataDDD(originalInd) = dataDDD(reorderedInd)
+    end if
+  end do
+
+
+  plotData = new( (/ nPlotYears, maxNddd /), typeof(data) )
+  plotDDD  = new( (/ nPlotYears, maxNddd /), typeof(dataDDD) )
+
+  do y = 0, dimsizes(i_plotYears)-1
+    isInYearInd := ind( dataYear.eq.i_plotYears(y) ) 
+    nCurrYear = dimsizes(isInYearInd)
+    plotDDD(y,1:nCurrYear) = dataDDD(isInYearInd)
+    plotData(y,1:nCurrYear) = (/ data(isInYearInd) /)
+
+    prevMaxInd = max( ind( dataYear.eq.i_plotYears(y)-1 ) )
+    if( .not.ismissing(prevMaxInd) ) then
+      plotData(y,0) = data(prevMaxInd)
+      plotDDD(y,0) = dataDDD(prevMaxInd)-365
+    end if
+
+    nextMinInd = min( ind( dataYear.eq.i_plotYears(y)+1 ) )
+    if( .not.ismissing(nextMinInd) ) then
+      plotData(y,nCurrYear+1) = data(nextMinInd)
+      plotDDD(y,nCurrYear+1) = dataDDD(nextMinInd)+365
+    end if
+  end do
+
+  firstMonth = toint( cd_string( i_startDDD, "%n" ) )
+  axisMonths = ispan( firstMonth, firstMonth + 12, 1 )
+  axisZeros  = axisMonths * 0
+  axisOnes   = axisZeros + 1
+  axisYears  = where( axisMonths.gt.12, 2, 1 )
+  axisMonths = where( axisMonths.gt.12, axisMonths-12, axisMonths )
+  axisTime = cd_inv_calendar( axisYears, axisMonths, axisOnes, \
+      axisZeros, axisZeros, axisZeros, i_startDDD@units, 0 )
+
+  ; Customize base plot
+  res = True
+  res@gsnDraw  = False
+  res@gsnFrame = False
+  res@gsnLeftString  = "" + data@long_name
+  res@gsnRightString = min(dataYear) + " - " + max(dataYear)
+  res@gsnRightString = min(dataYear) + " - " + 2014
+  res@tiYAxisString  = data@units
+  res@trXMinF = min(axisTime)
+  res@trXMaxF = max(axisTime)
+  res@tmXBLabelAngleF = 45
+  res@tmXBLabelJust = "CenterRight"
+  res@vpHeightF = 0.3
+
+  ; ...standardize the font sizes
+  res@gsnLeftStringFontHeightF  = i_fontHeightF
+  res@gsnRightStringFontHeightF = i_fontHeightF
+  res@tiMainFontHeightF         = i_fontHeightF * 1.3
+  res@tiYAxisFontHeightF        = i_fontHeightF
+  res@tiXAxisFontHeightF        = i_fontHeightF
+  res@tmXBLabelFontHeightF      = i_fontHeightF
+  res@tmYLLabelFontHeightF      = i_fontHeightF
+
+  ; ...make gridlines look nice
+  res@tmYMajorGridLineDashPattern = 2.
+  res@tmYMajorGridThicknessF      = 1.
+  res@tmYMajorGrid = True
+
+  ; ...set up an x-y line graph
+  res@xyMarkLineMode      = "Lines"
+  res@xyMonoDashPattern   = True
+  if( nPlotYears.le.6 ) then
+    res@pmLegendDisplayMode    = "Always"
+    res@xyLineColors     := (/ "red", "blue", "purple", "orange",  \
+        "magenta", "cyan", "seagreen", "brown" /)
+    res@xyLineThicknessF = 4
+  else
+    res@pmLegendDisplayMode    = "Never"
+    res@xyLineColor = "gray30"
+    res@xyLineThicknessF = 1
+  end if
+
+  ; ...set up the legend
+  res@lgPerimOn = True
+  res@lgPerimFill = 0
+  res@lgPerimFillColor = 0
+  res@lgPerimColor = 0
+  res@pmLegendSide           = "Top"
+  res@pmLegendZone           = 1
+  res@pmLegendParallelPosF   = 0.85
+  res@pmLegendOrthogonalPosF = 0.20
+  res@pmLegendWidthF         = 0.05
+  res@pmLegendHeightF        = 0.10
+  res@lgLabelFontHeightF     = i_fontHeightF * 0.7
+  res@xyExplicitLegendLabels = "   " + ( i_plotYears-1 ) + "/" \
+      + ( i_plotYears % 100 )
+  res@lgItemOrder = ispan( 0, nPlotYears-1, 1 )
+  res@lgItemOrder = res@lgItemOrder(::-1)
+
+  ; ...set up labels on the time axis
+  resTick = True
+  resTick@ttmFormat      = "%c"
+  resTick@ttmAxis        = "XB"
+  resTick@ttmMajorStride = 1
+  time_axis_labels( axisTime,  res, resTick )
+
+  retVal = gsn_csm_xy( io_wks, plotDDD, plotData, res )
+
+  return(retVal)
+end; draw_storage
+