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

grav = 9.80665
cp = 1004.6

dir = "./out"
dir = "./out-dwn1e1-lw"
dir = "./out-dwn1e1-0010-2000"
dir = "./out-dwn1e0-0010-2000"
dir = "./out-dwn1e0-0400-2000"
dirref = "."

sym    = "LW_Case25_Earth_Tropics_CO2-300ppmv"
symref = "Case25"
#sym    = "LW_Case09_Earth_MidLatSummer_CO2-300ppmvOnly"
#symref = "LBLRTM_Flux_Tendency_case09"
#sym    = "LW_Case19_Earth_MidLatSummer_H2OOnly"
#symref = "LBLRTM_Flux_Tendency_Case19"
#sym    = "LW_Case20_Earth_MidLatSummer_H2OOnly"
#symref = "LBLRTM_Flux_Tendency_Case20"
#sym    = "LW_Case23c_Earth_MidLatSummer_O3Only"
#symref = "LBLRTM_Flux_Tendency_Case23c"
#sym    = "LW_Case25_Earth_Tropics_CO2-300ppmv"
#symref = "LBLRTM_Flux_Tendency_Case25"
sym    = "LW_Case27_Earth_MidLatSummer_CO2-300ppmv"
symref = "LBLRTM_Flux_Tendency_Case27"


fn = sym + "_Flux.nc"
path = dir + "/" + fn
vname = 'RadUwFlux'
gpupfl = GPhys::NetCDF_IO.open( path, vname )
vname = 'RadDwFlux'
gpdnfl = GPhys::NetCDF_IO.open( path, vname )

plev = gpupfl.coord('Press').val
gpupfl = gpupfl.cut('Press'=>plev[0]..plev[-2])
gpdnfl = gpdnfl.cut('Press'=>plev[0]..plev[-2])

fn = sym + "_Tendency.nc"
path = dir + "/" + fn
vname = 'DTempDt'
gptend = GPhys::NetCDF_IO.open( path, vname )

plev = gptend.coord('Press').val
gptend = gptend.cut('Press'=>plev[0]..plev[-2])

fn = symref + ".nc"
dir = dirref
path = dir + "/" + fn
vname = 'RadUwFlux'
ref_gpupfl = GPhys::NetCDF_IO.open( path, vname )
vname = 'RadDwFlux'
ref_gpdnfl = GPhys::NetCDF_IO.open( path, vname )
#fn = symref + ".nc"
#vname = 'DTempDt'
#path = dir + "/" + fn
#ref_gptend = GPhys::NetCDF_IO.open( path, vname )
#ref_gptend = ref_gptend.copy

plev = ref_gpupfl.coord('Press').val
delp = plev[0..-2] - plev[1..-1]
ref_gptend = (ref_gpupfl-ref_gpdnfl)[0..-2] - (ref_gpupfl-ref_gpdnfl)[1..-1]
ref_gptend[true] = ref_gptend[true] / delp[true]
ref_gptend = ref_gptend * grav / cp


print "1: Display,  2: File\n"
citr = gets
citr = citr.chomp!
DCL.gropn(citr.to_i)
#DCL.gropn(1)

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

svx1 = 0.20; svx2 = 0.60; svy1 = 0.2; svy2 = 0.8

def graph_flux( svx1, svx2, svy1, svy2, gpupfl1, gpdnfl1, gpupfl2, gpdnfl2, flaglog, flagdiff )

  gpupfl1out = gpupfl1.copy
  gpdnfl1out = gpdnfl1.copy
  gpupfl2out = gpupfl2.copy
  gpdnfl2out = gpdnfl2.copy

  gpupfl1out.long_name = "flux"
  gpdnfl1out.long_name = "flux"

  x1 = 0; x2 = 500; y1 = 1.05e5; y2 = 0
#  x1 = 0; x2 = 1500; y1 = 1.05e5; y2 = 0
  itr = 1
  if flaglog then
    y2 = 1e1
    itr = 2
  end
  if flagdiff then
    x1 = -5; x2 = 5
#    x1 = -10; x2 = 10
#    x1 = -12; x2 = 12
    gpupfl2out.long_name = "flux diff."
    gpdnfl2out.long_name = "flux diff."

    gpupfl2out = gpupfl2out - gpupfl1out
    gpdnfl2out = gpdnfl2out - gpdnfl1out
    gpupfl1out = gpupfl1out - gpupfl1out
    gpdnfl1out = gpdnfl1out - gpdnfl1out
  end

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

  if !flagdiff then
    itype = 3
    gpout = gpupfl1out.cut('Press'=>y1..y2)
    GGraph.line( gpout, true , "exchange"=>true, "index"=>20, "type"=>itype, "title"=>"" )
    gpout = gpdnfl1out.cut('Press'=>y1..y2)
    GGraph.line( gpout, false, "exchange"=>true, "index"=>40, "type"=>itype, "title"=>"" )
  end

  if !flagdiff then
    flagfirst = false
  else
    flagfirst = true
  end
  itype = 1
  gpout = gpupfl2out.cut('Press'=>y1..y2)
  GGraph.line( gpout, flagfirst, "exchange"=>true, "index"=>20, "type"=>itype, "title"=>"" )
  gpout = gpdnfl2out.cut('Press'=>y1..y2)
  GGraph.line( gpout, false, "exchange"=>true, "index"=>40, "type"=>itype, "title"=>"" )

end

def graph_tendency( svx1, svx2, svy1, svy2, gp1, gp2, flaglog, flagdiff )

  gp1out = gp1.copy
  gp2out = gp2.copy

  if flagdiff then
    gp2out = gp2out - gp1out
    gp1out = gp1out - gp1out
  end

  gp1out = gp1out * 86400.0
  gp2out = gp2out * 86400.0

  gp1out.units = "K (day)-1"
  gp1out.long_name = "tendency"

  x1 = -5; x2 = 5; y1 = 1.05e5; y2 = 0
#  x1 = -20; x2 = 20; y1 = 1.05e5; y2 = 0
  itr = 1
  if flaglog then
    y2 = 1e1
    itr = 2
  end
  if flagdiff then
    x1 = -0.5; x2 = 0.5
#    x1 = -1  ; x2 = 1
#    x1 = -1.2; x2 = 1.2
#    x1 = -5  ; x2 = 5
    gp2out.long_name = "tendency diff."
  end

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

  if !flagdiff then
    itype = 3
    GGraph.line( gp1out, true , "exchange"=>true, "index"=>1, "title"=>"", "type"=>itype )
  end

  if !flagdiff then
    flagfirst = false
  else
    flagfirst = true
  end
  itype = 1
  GGraph.line( gp2out, flagfirst, "exchange"=>true, "index"=>1, "title"=>"", "type"=>itype )

end


# flux linear
graph_flux( svx1, svx2, svy1, svy2, ref_gpupfl, ref_gpdnfl, gpupfl, gpdnfl, false, false )

# flux diff. linear
graph_flux( svx1, svx2, svy1, svy2, ref_gpupfl, ref_gpdnfl, gpupfl, gpdnfl, false, true )

# flux log
graph_flux( svx1, svx2, svy1, svy2, ref_gpupfl, ref_gpdnfl, gpupfl, gpdnfl, true, false )

# flux diff. log
graph_flux( svx1, svx2, svy1, svy2, ref_gpupfl, ref_gpdnfl, gpupfl, gpdnfl, true, true )


# tendency linear
graph_tendency( svx1, svx2, svy1, svy2, ref_gptend, gptend, false, false )

# tendency diff. linear
graph_tendency( svx1, svx2, svy1, svy2, ref_gptend, gptend, false, true )

# tendency log
graph_tendency( svx1, svx2, svy1, svy2, ref_gptend, gptend, true, false )

# tendency diff. log
graph_tendency( svx1, svx2, svy1, svy2, ref_gptend, gptend, true, true )

DCL.grcls