#!/usr/bin/ruby
require "numru/ggraph"
include NumRu

# 120 bands
a_BandBnds = [10.0, 100.0, 200.0, 300.0, 400.0, 500.0, 600.0, 700.0, 800.0, 900.0, 1000.0, 1100.0, 1200.0, 1300.0, 1400.0, 1500.0, 1600.0, 1700.0, 1800.0, 1900.0, 2000.0, 2100.0, 2200.0, 2300.0, 2400.0, 2500.0, 2600.0, 2700.0, 2800.0, 2900.0, 3000.0, 3100.0, 3200.0, 3300.0, 3400.0, 3500.0, 3600.0, 3700.0, 3800.0, 3900.0, 4000.0, 4100.0, 4200.0, 4300.0, 4400.0, 4500.0, 4600.0, 4700.0, 4800.0, 4900.0, 5000.0, 5100.0, 5200.0, 5300.0, 5400.0, 5500.0, 5600.0, 5700.0, 5800.0, 5900.0, 6000.0, 6100.0, 6200.0, 6300.0, 6400.0, 6500.0, 6600.0, 6700.0, 6800.0, 6900.0, 7000.0, 7100.0, 7200.0, 7300.0, 7400.0, 7500.0, 7600.0, 7700.0, 8700.0, 9700.0, 10700.0, 11700.0, 12700.0, 13700.0, 14700.0, 15700.0, 16700.0, 17700.0, 18700.0, 19700.0, 20700.0, 21700.0, 22700.0, 23700.0, 24700.0, 25700.0, 26700.0, 27700.0, 28700.0, 29700.0, 30700.0, 31700.0, 32700.0, 33700.0, 34700.0, 35700.0, 36700.0, 37700.0, 38700.0, 39700.0, 40700.0, 41700.0, 42700.0, 43700.0, 44700.0, 45700.0, 46700.0, 47700.0, 48700.0, 49700.0, 50000.0]

grav           = 8.9     # gravitational acceleration (m -2)
mmw            = 44.0e-3 # mean molecular weight (kg mol-1)
dir            = "out"
ncfn_t         = "../prog01.0_mkprofile_ascii/out/Earth_ICRCCM_LW_Case27_MLS_CO2-300ppmv.nc"


# Below need not be changed.
nbands = a_BandBnds.size-1
a_iBands = [*(0..nbands)]
print nbands, " bands", "\n"
print a_BandBnds, "\n"

p_name    = "Press"
p_rms_min = 0e1
p_rms_max = 1e10

flag_diff = true

flag_disp_profile = true
flag_disp_profile = false


iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL.gropn(iws)

if ( flag_disp_profile ) then
  DCL.sldiv('y',4,2)
end
DCL.sgpset('lcntl',true)
DCL.sgpset('lfull',true)
#DCL.uzfact(1.5)
DCL.uzfact(0.6)            # font size
DCL.sgpset('lfprop',true)
DCL.sglset('lclip',true)


def graph_flux( svx1, svx2, svy1, svy2, x1, x2, gp, flaglog, label = "", flagxlog = false )

  if gp.instance_of?(Array) then
    gp_list = gp
  else
    gp_list = [gp]
  end


#  gpout = gp.copy
#
#  gpout.long_name = "flux"
#  gpout.units     = 'W m|-2"'


  y1 = 100.0e5; y2 = 0
#  y1 = 2.05e5; y2 = 0
  itr = 1
  if flaglog then
    y2 = 1e1
    y2 = 1e0
    itr = 2
    if flagxlog then
      itr = 4
    end
  else
    if flagxlog then
      itr = 3
    end
  end


  GGraph.set_fig 'itr'=>itr, 'viewport'=>[svx1,svx2,svy1,svy2], 'window'=>[x1,x2,y1,y2]

  flag_first = true
  gp_list.each_with_index { |gp_each, i|
    gpout = gp_each.cut('Press'=>y1..y2)
    GGraph.line( gpout, flag_first, "exchange"=>true, "index"=>(i+1)*10, "type"=>1, "title"=>"" )
    flag_first = false
  }


  rsize = DCL.sgqtxs()
  DCL.sgstxs( rsize*0.75 )
  DCL.sgstxi( 3 )
  DCL.sgstxc(-1)
  DCL.sglset('LCLIP', false)
  DCL.sgtxv( svx1, svy2*1.05, label )
  DCL.sglset('LCLIP', true)
  DCL.sgstxs( rsize )

end


def graph_tendency_log( svx1, svx2, svy1, svy2, x1, x2, gp, flaglog, label = "" )

  y1 = 100.0e5; y2 = 0
#  y1 = 2.05e5; y2 = 0
  itr = 1
  itr = 3
  if flaglog then
    y2 = 1e1
    y2 = 1e0
    itr = 2
    itr = 4
  end


  line_x_log( itr, svx1, svx2, svy1, svy2, x1, x2, y1, y2, gp )


  # panel symbol (a), (b), ...
  rsize = DCL.sgqtxs()
  DCL.sgstxs( rsize*0.75 )
  DCL.sgstxc(-1)
  DCL.sglset('LCLIP', false)
  DCL.sgtxv( svx1, svy2*1.05, label )
  DCL.sglset('LCLIP', true)
  DCL.sgstxs( rsize )

end


def line_x_log( itr, svx1, svx2, svy1, svy2, x1, x2, y1, y2, gp )


  if gp.instance_of?(Array) then
    gp_list = gp
  else
    gp_list = [gp]
  end

  gpn_list = []
  gpp_list = []
  gp_list.each { |gp_each|
    # negative
    gpn = gp_each.copy
    for k in 0..(gp_each.val.size-1)
      gpn[k] = [gp_each.val[k],-1e-10].min
    end
    # positive
    gpp = gp_each.copy
    for k in 0..(gp_each.val.size-1)
      gpp[k] = [gp_each.val[k],1e-10].max
    end
    gpn_list << gpn
    gpp_list << gpp
  }


  DCL.sgfrm

  svx1_local = svx1
  svx2_local = svx1 + ( svx2 - svx1 ) / 2 - ( svx2 - svx1 ) / 2 * 0.02
  x1_local   = -x2
  x2_local   = -x1
  #
  DCL.grfig
  #
  DCL.sgsvpt(svx1_local,svx2_local,svy1,svy2)
  DCL.sgswnd(x1_local,x2_local,y1,y2)
  DCL.grstrn(itr)
  DCL.grstrf
  #
  lidx = DCL.uuqlni()
  gpn_list.each_with_index { |gpn, i|
    DCL.uuslni((i+1)*10)
    DCL.uulin(gpn.val, gpn.coord('Press').val)
  }
  DCL.uuslni(lidx)
  #
#  DCL.ussttl( gp2outsign.long_name.to_s, '', gp2outsign.coord('Press').long_name, gp2outsign.coord('Press').units.to_s )
  #
  DCL.uzlset( 'labelxt', false )
  DCL.uzlset( 'labelxb', false )
  DCL.uzlset( 'labelyl', true  )
  DCL.uzlset( 'labelyr', false )
#  DCL.uscset( 'cxside', 'bt' )
#  DCL.uscset( 'cyside', 'lr' )
  DCL.usdaxs
  # numeric labels for x
  flabels=[]
  clabels=[]
  log10(x2).to_i.step( log10(x1).to_i, -4 ) do |i|
    flabels << -10**i
    clabels << '-10|' + i.to_s + '"'
  end
  DCL.uxplbl('b',1,flabels,clabels,4)


  svx1_local = svx1 + ( svx2 - svx1 ) / 2 + ( svx2 - svx1 ) / 2 * 0.02
  svx2_local = svx2
  x1_local   = x1
  x2_local   = x2
  #
  DCL.grfig
  #
  DCL.sgsvpt(svx1_local,svx2_local,svy1,svy2)
  DCL.sgswnd(x1_local,x2_local,y1,y2)
  DCL.grstrn(itr)
  DCL.grstrf
  #
  lidx = DCL.uuqlni()
  gpp_list.each_with_index { |gpp, i|
    DCL.uuslni((i+1)*10)
    DCL.uulin(gpp.val, gpp.coord('Press').val)
  }
  DCL.uuslni(lidx)
  #
#  DCL.ussttl( '', gp2outsign.units.to_s, gp2outsign.coord('Press').long_name, gp2outsign.coord('Press').units.to_s )
  DCL.uzlset( 'labelxt', false )
  DCL.uzlset( 'labelxb', false )
  DCL.uzlset( 'labelyl', false )
  DCL.uzlset( 'labelyr', false )
  DCL.usdaxs
  # numeric labels for x
  flabels=[]
  clabels=[]
  log10(x2).to_i.step( log10(x1).to_i, -4 ) do |i|
    flabels << 10**i
    clabels << '10|' + i.to_s + '"'
  end
  DCL.uxplbl('b',1,flabels,clabels,4)
  #
  # Label
  DCL.sgsvpt(svx1,svx2,svy1,svy2)
  DCL.grstrf
  rsize = DCL.sgqtxs()
  DCL.sgstxs( rsize*0.5 )
  rsizec = DCL.uzpget( 'rsizec1' )
  DCL.uzpset( 'rsizec1', rsizec * 0.75 )
  DCL.uxsttl('b',"("+gpp_list[0].units.to_s+")",1)
  DCL.uzpset( 'rsizec1', rsizec )
  DCL.sgstxs( rsize )
  DCL.uxsttl('b',gpp_list[0].long_name.to_s,0)
  DCL.uysttl('l',' ',0)
  DCL.uysttl('l',gpp_list[0].coord('Press').long_name.to_s,0)
#  DCL.ussttl( '', gp2outsign.units.to_s, gp2outsign.coord('Press').long_name, gp2outsign.coord('Press').units.to_s )

  # reset setting
  DCL.uzlset( 'labelxb', true )
  DCL.uzlset( 'labelyl', true )

end


# calculate height
gasconst = 8.314 / mmw
gp_temp  = GPhys::NetCDF_IO.open( ncfn_t, "Temp" )
na_temp  = gp_temp.val
na_press = gp_temp.coord('Press').val
if ( na_press[-1] == 0.0 ) then
  na_press[-1] = na_press[-2] / 2
end
# dp/dz = -rho g
# dz = - 1/(rho g) dp = - p / (rho g) d(ln(p))
#    = - rho R T / (rho g) d(ln(p)) = - R T / g d(ln(p))
gpz = gp_temp.copy
gpz.units = "m"
kmax = gp_temp.val.size
k = 0
gpz[k] = 0.0
for k in 1..(kmax-1)
  gpz[k] = gpz.val[k-1] -
            gasconst * ( na_temp[k] + na_temp[k-1] ) / 2 / grav *
            log( na_press[k] / na_press[k-1] )
end



pr_upflx  = []
pr_dnflx  = []
pr_flxcnv = []
pr_tend   = []
sr_upflx  = []
sr_dnflx  = []
sr_flxcnv = []
sr_tend   = []

a_iBands.each{ |iBand|

  ncfn_flux         = dir+"/"+"Flux-"    +iBand.to_s.rjust(3,"0")+"-10.nc"
  ncfn_tendency     = dir+"/"+"Tendency-"+iBand.to_s.rjust(3,"0")+"-10.nc"
  ncfn_flux_ref     = dir+"/"+"Flux-"    +iBand.to_s.rjust(3,"0")+"-LBL.nc"
  ncfn_tendency_ref = dir+"/"+"Tendency-"+iBand.to_s.rjust(3,"0")+"-LBL.nc"

  icharnum = "a".ord - 1
  #
  exts = ['LW', 'SW']
  for iext in 0..(exts.size-1)
    ext = exts[iext]

#    if ( a_BandBnds.size > 2 ) then
#      if ( iBand == 0 ) then
#        print ext, ", ", iBand, ", ", a_BandBnds[0], "-", a_BandBnds[-1], "\n"
#      else
#        print ext, ", ", iBand, ", ", a_BandBnds[iBand-1], "-", a_BandBnds[iBand], "\n"
#      end
#    else
#      print ext, ", ", iBand, "\n"
#    end

    # k-dist
    #   UwFlux
    path = ncfn_flux
    vname = 'RadUwFlux' + ext
    gpupfl = GPhys::NetCDF_IO.open( path, vname )
    gp = gpupfl.copy
    z = gp.axis("Press").pos ; z.long_name = 'pressure' ; gp.axis("Press").set_pos(z)
    gpupfl = gp
    #   DwFlux
    path = ncfn_flux
    vname = 'RadDwFlux' + ext
    gpdnfl = GPhys::NetCDF_IO.open( path, vname )
    gp = gpdnfl.copy
    z = gp.axis("Press").pos ; z.long_name = 'pressure' ; gp.axis("Press").set_pos(z)
    gpdnfl = gp
    #   FluxConv
    kmax = gpupfl.coord('Press').val.size
    gpnetfl = gpupfl - gpdnfl
    gpnetflconv = - ( gpnetfl[1..(kmax-1)] - gpnetfl[0..(kmax-2)] ) / ( gpz[1..(kmax-1)] - gpz[0..(kmax-2)] )
    #   Tendency
    path = ncfn_tendency
    vname = 'DTempDt' + ext
    gpdtdt = GPhys::NetCDF_IO.open( path, vname )
    gp = gpdtdt.copy
    z = gp.axis("Press").pos ; z.long_name = 'pressure' ; gp.axis("Press").set_pos(z)
    gpdtdt = gp

    # LBL
    #   UwFlux
    path = ncfn_flux_ref
    vname = 'RadUwFlux' + ext
    ref_gpupfl = GPhys::NetCDF_IO.open( path, vname )
    #   DwFlux
    path = ncfn_flux_ref
    vname = 'RadDwFlux' + ext
    ref_gpdnfl = GPhys::NetCDF_IO.open( path, vname )
    #   FluxConv
    kmax = ref_gpupfl.coord('Press').val.size
    ref_gpnetfl = ref_gpupfl - ref_gpdnfl
    ref_gpnetflconv = - ( ref_gpnetfl[1..(kmax-1)] - ref_gpnetfl[0..(kmax-2)] ) / ( gpz[1..(kmax-1)] - gpz[0..(kmax-2)] )
    #   Tendency
    path = ncfn_tendency_ref
    vname = 'DTempDt' + ext
    ref_gpdtdt = GPhys::NetCDF_IO.open( path, vname )
    ref_gpdtdt = ref_gpdtdt.copy

    svx1 = 0.175; svx2 = 0.65; svy1 = 0.2; svy2 = 0.8
    svx1 = 0.175; svx2 = 0.63; svy1 = 0.2; svy2 = 0.8

    # Uw flux

    # flux diff.
    if flag_diff then
      gpupfl = gpupfl - ref_gpupfl
      gpdnfl = gpdnfl - ref_gpdnfl
      gpupfl.long_name = "up.flux diff."
      gpdnfl.long_name = "dn.flux diff."
    else
      gpupfl.long_name = "up.flux"
      gpdnfl.long_name = "dn.flux"
    end
    gpupfl.units = 'W m|-2"'
    gpdnfl.units = 'W m|-2"'
    if flag_disp_profile then
      if flag_diff then
        # difference
        x1 = -4; x2 = 4
      else
        if exts[iext] == 'LW' then
          x1 = 0; x2 = 2e4
        else
          x1 = 0; x2 = 1e3
        end
      end
      #
      icharnum += 1
      label = "("+icharnum.chr+")"
      label += " " + ext
      if flag_diff then
        gp_list = [gpupfl]
      else
        gp_list = [ref_gpupfl,gpupfl]
      end
      graph_flux( svx1, svx2, svy1, svy2, x1, x2, gp_list, true, label )
    end
    #
    # Print RMS error
    na_values = gpupfl.cut(p_name=>p_rms_min..p_rms_max).val
    rms = sqrt( ( na_values**2 ).mean )
#    print 'upfl_rms    ', rms, "\n"
    if ( ext == 'LW' ) then
      pr_upflx  << rms
    else
      sr_upflx  << rms
    end

    #
    # Dw flux
    if flag_disp_profile then
      if flag_diff then
        # difference
        x1 = -4; x2 = 4
      else
        if exts[iext] == 'LW' then
          x1 = 0; x2 = 2e4
        else
          x1 = 0; x2 = 1e3
        end
      end
      icharnum += 1
      label = "("+icharnum.chr+")"
      label += " " + ext
      if flag_diff then
        gp_list = [gpdnfl]
      else
        gp_list = [ref_gpdnfl,gpdnfl]
      end
      graph_flux( svx1, svx2, svy1, svy2, x1, x2, gp_list, true, label )
    end
    #
    # Print RMS error
    na_values = gpdnfl.cut(p_name=>p_rms_min..p_rms_max).val
    rms = sqrt( ( na_values**2 ).mean )
#    print 'dnfl_rms    ', rms, "\n"
    if ( ext == 'LW' ) then
      pr_dnflx  << rms
    else
      sr_dnflx  << rms
    end


    # Flux convergence

    # tendency diff. log
    if flag_diff then
      gpnetflconv = gpnetflconv - ref_gpnetflconv
      gpnetflconv.long_name = "flux conv. diff."
    else
      gpnetflconv.long_name = "flux conv."
    end
    gpnetflconv.units = 'W m|-3"'
    #
    if flag_disp_profile then
      if flag_diff then
        # difference
        x1 = -1e-3; x2 = 1e-3
      else
        x1 = -1.5e-2; x2 = 1.5e-2
      end
      #
      x2 = gpnetflconv.val.abs.max*1.05
      x2 = 4e-4
      x1 = -x2
      #
      icharnum += 1
      label = "("+icharnum.chr+")"
      label += " " + ext
      if flag_diff then
        gp_list = [gpnetflconv]
      else
        gp_list = [ref_gpnetflconv,gpnetflconv]
      end
#    graph_tendency_log( svx1, svx2, svy1, svy2, x1, x2, gptend_list, true, label )
      graph_flux( svx1, svx2, svy1, svy2, x1, x2, gp_list, true, label )
    end
    #
    # Print RMS error
    na_values = gpnetflconv.cut(p_name=>p_rms_min..p_rms_max).val
    rms = sqrt( ( na_values**2 ).mean )
#    print 'flxcnv_rms  ', rms, "\n"
    if ( ext == 'LW' ) then
      pr_flxcnv << rms
    else
      sr_flxcnv << rms
    end


    # Tendency

    # tendency diff. log
    if flag_diff then
      gpdtdt = gpdtdt - ref_gpdtdt
      gpdtdt.long_name = "tendency diff."
    else
      gpdtdt.long_name = "tendency"
    end
    gpdtdt.units = 'K s|-1"'
    #
    if flag_disp_profile then
      if flag_diff then
        # difference
        x1 = -1e-4; x2 = 1e-4
      else
        x1 = -2e-4; x2 = 2e-4
      end
      #
      x2 = gpdtdt.val.abs.max*1.05
      x1 = -x2
      #
      icharnum += 1
      label = "("+icharnum.chr+")"
      label += " " + ext
      if flag_diff then
        gp_list = [gpdtdt]
      else
        gp_list = [ref_gpdtdt,gpdtdt]
      end
#    graph_tendency_log( svx1, svx2, svy1, svy2, x1, x2, gptend_list, true, label )
      graph_flux( svx1, svx2, svy1, svy2, x1, x2, gp_list, true, label )
    end
    #
    # Print RMS error
    na_values = gpdtdt.cut(p_name=>p_rms_min..p_rms_max).val
    rms = sqrt( ( na_values**2 ).mean )
#    print 'tend_rms    ', rms, "\n"
    if ( ext == 'LW' ) then
      pr_tend   << rms
    else
      sr_tend   << rms
    end

#  DCL.grfrm

  end

}


a_iBands.each{ |iBand|
  if ( a_BandBnds.size > 2 ) then
    if ( iBand == 0 ) then
      print iBand, ", ", a_BandBnds[0], "-", a_BandBnds[-1], "\n"
    else
      print iBand, ", ", a_BandBnds[iBand-1], "-", a_BandBnds[iBand], "\n"
    end
  else
    print iBand, "\n"
  end
  print "           Planetary Rad.      Solar Rad.", "\n"
  str_upflx  = sprintf("%10.5e    |    %10.5e", pr_upflx [iBand], sr_upflx [iBand])
  str_dnflx  = sprintf("%10.5e    |    %10.5e", pr_dnflx [iBand], sr_dnflx [iBand])
  str_flxcnv = sprintf("%10.5e    |    %10.5e", pr_flxcnv[iBand], sr_flxcnv[iBand])
  str_tend   = sprintf("%10.5e    |    %10.5e", pr_tend  [iBand], sr_tend  [iBand])
  print '  ', 'upflx  : ', str_upflx , "\n"
  print '  ', 'dnflx  : ', str_dnflx , "\n"
  print '  ', 'flxcnv : ', str_flxcnv, "\n"
  print '  ', 'tend   : ', str_tend  , "\n"
}

a_X = []
for i in 0..(a_BandBnds.size-1)
  a_X << a_BandBnds[i] ; a_X << a_BandBnds[i]
end
a_rmse = pr_upflx
a_Y = []
for i in 0..(a_rmse.size-1)
  a_Y << a_rmse[i] ; a_Y << a_rmse[i]
end


def graph_rmse_spe( vname, units, a_iBands, a_BandBnds, a_RMSE_list )

  maxval = 0.0
  for i in 0..(a_RMSE_list.size-1)
    a_RMSE = a_RMSE_list[i]
    maxval = [maxval, a_RMSE.max].max
  end
  flag_first = true
  for i in 0..(a_RMSE_list.size-1)
    a_RMSE = a_RMSE_list[i]
    graph_rmse_spe_core( vname, units, flag_first, a_iBands, a_BandBnds, a_RMSE, true, maxval )
    flag_first = false
  end

end

def graph_rmse_spe_core( vname, units, flag_first, a_iBands, a_BandBnds, a_RMSE, flag_bandnum = false, maxval )

  itr = 3

  x1 = a_BandBnds[0] * 0.9
  x2 = a_BandBnds[-1] * 1.1
  #y1 = gp.val.min
  y1 = 0
#  y2 = gp.val.max * 1.1
  y2 = maxval * 1.1
  GGraph.set_fig 'itr'=>itr, 'viewport'=>[0.1,0.9,0.1,0.65], 'window'=>[x1,x2,y1,y2]


  if ( flag_first ) then
    flag_first_local = true
  else
    flag_first_local = false
  end
  a_iBands.each{ |iBand|
    a_X = []
    a_Data = []
    if ( iBand == 0 ) then
      a_X    << a_BandBnds[0]
      a_X    << a_BandBnds[-1]
      a_Data << a_RMSE[0]
      a_Data << a_RMSE[0]
    else
      a_X    << a_BandBnds[iBand-1]
      a_X    << a_BandBnds[iBand-1]
      a_X    << a_BandBnds[iBand]
      a_X    << a_BandBnds[iBand]
      a_Data << 0.0
      a_Data << a_RMSE[iBand]
      a_Data << a_RMSE[iBand]
      a_Data << 0.0
    end
    va_X = VArray.new( NArray.to_na(a_X),
                       {"long_name"=>"wavenumber", "units"=>"cm-1"},
                       "wn" )
    ax_X = Axis.new.set_pos(va_X)
    va_Data = VArray.new( NArray.to_na(a_Data), 
                          {"long_name"=>"RMSE"+" of "+vname, "units"=>units},
                          "RMSE" )
    gp = GPhys.new( Grid.new(ax_X), va_Data )
    if ( flag_first ) then
      index = 20
    else
      index = 40
    end
    if ( iBand == 0 ) then
      itype = 2
    else
      itype = 1
    end
    GGraph.line gp, flag_first_local, 'exchange'=>false, 'index'=>index, 'type'=>itype
    flag_first_local = false

    # band number label
    if ( flag_bandnum ) then
      iindex = DCL.sgqtxi()
      DCL.sgstxi(index)
      rsize = DCL.sgqtxs()
      DCL.sgstxs( rsize*0.25 )
      DCL.sgstxc(0.5)
      if ( iBand == 0 ) then
        wns = a_BandBnds[0] ; wne = a_BandBnds[-1]
        str = "total"
      else
        wns = a_BandBnds[iBand-1] ; wne = a_BandBnds[iBand]
        str = iBand.to_s
      end
      if ( itr == 3 ) then
        x = sqrt( wns * wne )
      else
        x = ( wns + wne ) / 2
      end
      y = a_RMSE[iBand] + (y2-y1)*0.025
      DCL.sgtxu( x, y, str )
      DCL.sgstxc(0)
      DCL.sgstxs( rsize )
      DCL.sgstxi(iindex)
    end
  }

end

unless ( flag_disp_profile ) then
  a_rmse_list = [pr_upflx, sr_upflx]
  graph_rmse_spe( 'UpFlux', 'W m-2', a_iBands, a_BandBnds, a_rmse_list )
  a_rmse_list = [pr_dnflx, sr_dnflx]
  graph_rmse_spe( 'DnFlux', 'W m-2', a_iBands, a_BandBnds, a_rmse_list )
  a_rmse_list = [pr_flxcnv, sr_flxcnv]
  graph_rmse_spe( 'FluxConv', 'W m-3', a_iBands, a_BandBnds, a_rmse_list )
  a_rmse_list = [pr_tend, sr_tend]
  graph_rmse_spe( 'Tendency', 'K s-1', a_iBands, a_BandBnds, a_rmse_list )
end

DCL.grcls