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


grav           = 9.8     # gravitational acceleration (m -2)
mmw            = 29.0e-3 # mean molecular weight (kg mol-1)
ncfn_t         = "lbl_k-dist_test_TRAPPIST-1/prog01.0_mkprofile_ascii/out/Earth_ICRCCM_LW_Case27_MLS_CO2-300ppmv.nc"
ncfn1_flux     = "lbl_k-dist_test_TRAPPIST-1/prog03.0_calc_rte/out/TRAPPIST-1_planet_Flux.nc"
ncfn1_tendency = "lbl_k-dist_test_TRAPPIST-1/prog03.0_calc_rte/out/TRAPPIST-1_planet_Tendency.nc"
ncfn2_flux     = "lbl_k-dist_test_TRAPPIST-1_dwn1e2/prog03.0_calc_rte/out/TRAPPIST-1_planet_Flux.nc"
ncfn2_tendency = "lbl_k-dist_test_TRAPPIST-1_dwn1e2/prog03.0_calc_rte/out/TRAPPIST-1_planet_Tendency.nc"


# 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


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

DCL.sldiv('y',4,2)
DCL.sgpset('lcntl',true)
DCL.sgpset('lfull',true)
DCL.uzfact(1.5)
DCL.sgpset('lfprop',true)
DCL.sglset('lclip',true)


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

  y1 = 1.1e5; 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
  gpout = gp.cut('Press'=>y1..y2)
  GGraph.line( gpout, flag_first, "exchange"=>true, "index"=>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


icharnum = "a".ord - 1
#
exts = ['LW', 'SW']
for iext in 0..(exts.size-1)
  ext = exts[iext]
  if ( iext == 0 ) then
    label_ext = 'PR'
  else
    label_ext = 'SR'
  end

  gpupfl_list      = []
  gpdnfl_list      = []
  gpnetflconv_list = []
  gpdtdt_list      = []

  i = 0
  ncfn_flux     = ncfn1_flux
  ncfn_tendency = ncfn1_tendency

  # UwFlux, DwFlux

  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
  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

  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


  ncfn_flux_ref     = ncfn2_flux
  ncfn_tendency_ref = ncfn2_tendency
  #
  path = ncfn_flux_ref
  vname = 'RadUwFlux' + ext
  ref_gpupfl = GPhys::NetCDF_IO.open( path, vname )
  vname = 'RadDwFlux' + ext
  ref_gpdnfl = GPhys::NetCDF_IO.open( path, vname )
  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)] )
  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.63; svy1 = 0.2; svy2 = 0.8

  # Uw flux

  # difference
  x1 = -4; x2 = 4
#      x1 = -50; x2 = 50

  # flux diff.
  gpupfl = gpupfl - ref_gpupfl
  gpdnfl = gpdnfl - ref_gpdnfl
  gpupfl.long_name = "up.flux diff."
  gpdnfl.long_name = "dn.flux diff."
  gpupfl.units = 'W m|-2"'
  gpdnfl.units = 'W m|-2"'

  icharnum += 1
  label = "("+icharnum.chr+")"
  label += " " + label_ext
  mkgraph( svx1, svx2, svy1, svy2, x1, x2, gpupfl, true, label )
  #
  # Print RMS error
#  print "-----\n"
#  print "# ", Dir.pwd, "\n"
#  print "# This is calculated by disp_comp_fluxfluxconv_multi_1panel_v5.rb.\n"
#  print "# ", label, "\n"
#  print "# file     upfl_rms\n"
  na_upfl = gpupfl.val
  na_upfl = na_upfl**2
  rmse    = na_upfl.mean
  rmse    = sqrt( rmse )

  if ( ext == 'LW' ) then
    pr_upflx = rmse
  else
    sr_upflx = rmse
  end
  #
  # difference
  x1 = -4; x2 = 4
  #
  # Dw flux

  icharnum += 1
  label = "("+icharnum.chr+")"
  label += " " + label_ext
  if exts[iext] == 'LW' then
    mkgraph( svx1, svx2, svy1, svy2, x1, x2, gpdnfl, true, label )
  else
    mkgraph( svx1, svx2, svy1, svy2, x1, x2, gpdnfl, true, label )
  end
  #
  # Print RMS error
#  print "-----\n"
#  print "# ", Dir.pwd, "\n"
#  print "# This is calculated by disp_comp_fluxfluxconv_multi_1panel_v5.rb.\n"
#  print "# ", label, "\n"
#  print "# file     dnfl_rms\n"
  na_dnfl = gpdnfl.val
  na_dnfl  = na_dnfl**2
  rmse = na_dnfl.mean
  rmse = sqrt( rmse )

  if ( ext == 'LW' ) then
    pr_dnflx = rmse
  else
    sr_dnflx = rmse
  end


  # Flux convergence

  icharnum += 1
  # tendency diff. log
  gpnetflconv = gpnetflconv - ref_gpnetflconv
  gpnetflconv.long_name = "flux conv. diff."
  gpnetflconv.units = 'W m|-3"'
  # difference
  x1 = -1e-3; x2 = 1e-3
  label = "("+icharnum.chr+")"
  label += " " + label_ext
  mkgraph( svx1, svx2, svy1, svy2, x1, x2, gpnetflconv, true, label )
  #
  # Print RMS error
#  print "-----\n"
#  print "# ", Dir.pwd, "\n"
#  print "# This is calculated by disp_comp_fluxfluxconv_multi_1panel_v5.rb.\n"
#  print "# ", label, "\n"
#  print "# file     flxcnv_rms\n"
  na_tend = gpnetflconv.val
  na_tend  = na_tend**2
  rmse = na_tend.mean
  rmse = sqrt( rmse )

  if ( ext == 'LW' ) then
    pr_flxcnv = rmse
  else
    sr_flxcnv = rmse
  end


  # Tendency

  icharnum += 1
  # tendency diff. log
  gpdtdt = gpdtdt - ref_gpdtdt
  gpdtdt.long_name = "tendency diff."
  gpdtdt.units = 'K s|-1"'
  # difference
  x1 = -1e-4; x2 = 1e-4
  #
  label = "("+icharnum.chr+")"
  label += " " + label_ext
  mkgraph( svx1, svx2, svy1, svy2, x1, x2, gpdtdt, true, label )
  #
  # Print RMS error
#  print "-----\n"
#  print "# ", Dir.pwd, "\n"
#  print "# This is calculated by disp_comp_fluxfluxconv_multi_1panel_v5.rb.\n"
#  print "# ", label, "\n"
#  print "# file     tend_rms\n"
  na_tend = gpdtdt.val
  na_tend  = na_tend**2
  rmse = na_tend.mean
  rmse = sqrt( rmse )

  if ( ext == 'LW' ) then
    pr_tend = rmse
  else
    sr_tend = rmse
  end

end

DCL.grcls

print "RMSE summary", "\n"
print " 1st file  : ", ncfn1_flux, "\n"
print "           : ", ncfn1_tendency, "\n"
print " 2nd file  : ", ncfn2_flux, "\n"
print "           : ", ncfn2_tendency, "\n"
print " variable : Planetary Rad.      Solar Rad.", "\n"
str_upflx  = sprintf("%10.5e    |    %10.5e W m-2", pr_upflx , sr_upflx )
str_dnflx  = sprintf("%10.5e    |    %10.5e W m-2", pr_dnflx , sr_dnflx )
str_flxcnv = sprintf("%10.5e    |    %10.5e W m-3", pr_flxcnv, sr_flxcnv)
str_tend   = sprintf("%10.5e    |    %10.5e K s-1", pr_tend  , sr_tend  )
print '  ', 'upflx   : ', str_upflx , "\n"
print '  ', 'dnflx   : ', str_dnflx , "\n"
print '  ', 'flxcnv  : ', str_flxcnv, "\n"
print '  ', 'tend    : ', str_tend  , "\n"