# -*- coding: euc-jp -*-
# Title: Ruby script drawing contour map for deepconv/arare5 output data 
#
# History: 2011/09/27 (Masatsugu Odaka)
#
require "numru/ggraph"

module NumRu
  class GPhys
    def deriv(dim)
      dim = dim_index(dim)    # input dim can be a String or an Integer
      x = axis(dim).to_gphys
      a = [true]*dim + [1..-1,false]
      b = [true]*dim + [0..-2,false]
      dydx = ( self[*a] - self[*b] ) / ( x[1..-1] - x[0..-2] )
      xi = ( x[1..-1] + x[0..-2] ) / 2
      dydx.axis(dim).set_pos(xi.coord(0))
      dydx.long_name = "d #{name}/d #{x.name}"
      dydx
    end
  end
end


include NumRu

filedir = "time_000180000-000280000"
filehead = "BS1998_"
cuttime = 180000..280000
puttime = 50000
zcut = 10000
kappa = 155.0
comax =5.0e-7
comin = -5.0e-7
ins = 2
#comax =0.000010
#comin = -0.0000030
#file0 = 'BS1998_restart_rank000000.nc'
dens_file = 'restart_long-DensBZ.nc'
ptemp0_file = "../#{filedir}/#{filehead}PTemp_rank000000.nc"
ptemp1_file = "../#{filedir}/#{filehead}PTemp_rank000001.nc"
ptemp2_file = "../#{filedir}/#{filehead}PTemp_rank000002.nc"
ptemp3_file = "../#{filedir}/#{filehead}PTemp_rank000003.nc"
ptemp4_file = "../#{filedir}/#{filehead}PTemp_rank000004.nc"
ptemp5_file = "../#{filedir}/#{filehead}PTemp_rank000005.nc"

adv0_file = "../#{filedir}/#{filehead}PTempAdv_rank000000.nc"
adv1_file = "../#{filedir}/#{filehead}PTempAdv_rank000001.nc"
adv2_file = "../#{filedir}/#{filehead}PTempAdv_rank000002.nc"
adv3_file = "../#{filedir}/#{filehead}PTempAdv_rank000003.nc"
adv4_file = "../#{filedir}/#{filehead}PTempAdv_rank000004.nc"
adv5_file = "../#{filedir}/#{filehead}PTempAdv_rank000005.nc"

turb0_file = "../#{filedir}/#{filehead}PTempTurb_rank000000.nc"
turb1_file = "../#{filedir}/#{filehead}PTempTurb_rank000001.nc"
turb2_file = "../#{filedir}/#{filehead}PTempTurb_rank000002.nc"
turb3_file = "../#{filedir}/#{filehead}PTempTurb_rank000003.nc"
turb4_file = "../#{filedir}/#{filehead}PTempTurb_rank000004.nc"
turb5_file = "../#{filedir}/#{filehead}PTempTurb_rank000005.nc"

sfc0_file = "../#{filedir}/#{filehead}PTempSfc_rank000000.nc"
sfc1_file = "../#{filedir}/#{filehead}PTempSfc_rank000001.nc"
sfc2_file = "../#{filedir}/#{filehead}PTempSfc_rank000002.nc"
sfc3_file = "../#{filedir}/#{filehead}PTempSfc_rank000003.nc"
sfc4_file = "../#{filedir}/#{filehead}PTempSfc_rank000004.nc"
sfc5_file = "../#{filedir}/#{filehead}PTempSfc_rank000005.nc"

rad0_file = "../#{filedir}/#{filehead}PTempRad_rank000000.nc"
rad1_file = "../#{filedir}/#{filehead}PTempRad_rank000001.nc"
rad2_file = "../#{filedir}/#{filehead}PTempRad_rank000002.nc"
rad3_file = "../#{filedir}/#{filehead}PTempRad_rank000003.nc"
rad4_file = "../#{filedir}/#{filehead}PTempRad_rank000004.nc"
rad5_file = "../#{filedir}/#{filehead}PTempRad_rank000005.nc"

velx0_file = "../#{filedir}/#{filehead}VelX_rank000000.nc"
velx1_file = "../#{filedir}/#{filehead}VelX_rank000001.nc"
velx2_file = "../#{filedir}/#{filehead}VelX_rank000002.nc"
velx3_file = "../#{filedir}/#{filehead}VelX_rank000003.nc"
velx4_file = "../#{filedir}/#{filehead}VelX_rank000004.nc"
velx5_file = "../#{filedir}/#{filehead}VelX_rank000005.nc"

velz0_file = "../#{filedir}/#{filehead}VelZ_rank000000.nc"
velz1_file = "../#{filedir}/#{filehead}VelZ_rank000001.nc"
velz2_file = "../#{filedir}/#{filehead}VelZ_rank000002.nc"
velz3_file = "../#{filedir}/#{filehead}VelZ_rank000003.nc"
velz4_file = "../#{filedir}/#{filehead}VelZ_rank000004.nc"
velz5_file = "../#{filedir}/#{filehead}VelZ_rank000005.nc"

ptemp  = GPhys::IO.open([ptemp0_file, \
                         ptemp1_file, \
                         ptemp2_file, \
                         ptemp3_file, \
                         ptemp4_file, \
                         ptemp5_file ], \
                        'PTemp')

adv    = GPhys::IO.open([adv0_file, \
                         adv1_file, \
                         adv2_file, \
                         adv3_file, \
                         adv4_file, \
                         adv5_file ], \
                        'PTempAdv')


turb   = GPhys::IO.open([turb0_file, \
                         turb1_file, \
                         turb2_file, \
                         turb3_file, \
                         turb4_file, \
                         turb5_file ], \
                        'PTempTurb')

rad    = GPhys::IO.open([rad0_file, \
                         rad1_file, \
                         rad2_file, \
                         rad3_file, \
                         rad4_file, \
                         rad5_file ], \
                        'PTempRad')

sfc    = GPhys::IO.open([sfc0_file, \
                         sfc1_file, \
                         sfc2_file, \
                         sfc3_file, \
                         sfc4_file, \
                         sfc5_file ], \
                        'PTempSfc')

velx   = GPhys::IO.open([velx0_file, \
                         velx1_file, \
                         velx2_file, \
                         velx3_file, \
                         velx4_file, \
                         velx5_file ], \
                        'VelX')
velz   = GPhys::IO.open([velz0_file, \
                         velz1_file, \
                         velz2_file, \
                         velz3_file, \
                         velz4_file, \
                         velz5_file ], \
                        'VelZ')

densBZ  = GPhys::IO.open(dens_file, 'DensBZ')

turb = turb + sfc


# nonconv นเ

  rho_u = densBZ * velx
  rho_w = densBZ * velz
  
  dx_rho_u = rho_u.deriv('x')
  dz_rho_w = rho_w.deriv('z')
  
  conv = ptemp[1..1199,0..-1,1..99,0..-1] * ( dx_rho_u[0..-1,0..-1,1..99,0..-1] +  dz_rho_w[1..1199,0..-1,0..-1,0..-1] )
  
  conv_XYZmean = conv.mean('x').mean('y').mean('z')
  conv_XYmean = conv.mean('x').mean('y')
  conv_XYmeanZcut = conv_XYmean.cut('z'=>zcut)

  p conv_XYZmean.cut('t'=>0.0).val.to_f
  p conv_XYZmean.cut('t'=>80000).val.to_f


# rho theta

  rho_theta = ptemp * densBZ
  #rho_theta = ptemp[1..1199,0..-1,1..99,0..-1] * densBZ[1..1199,0..-1,1..99]
  dt_rho_theta = rho_theta.deriv('t')
  dt_rho_theta_XYZmean = dt_rho_theta.mean('x').mean('y').mean('z')
  dt_rho_theta_XYmean = dt_rho_theta.mean('x').mean('y')
  dt_rho_theta_XYmeanZcut = dt_rho_theta_XYmean.cut('z'=>zcut)

# turb
  dxdx_theta = kappa * ptemp.deriv('x').deriv('x') * densBZ[2..1199,0..-1,0..-1]
  dz_theta = ptemp.deriv('z') * densBZ[0..-1,0..-1,1..99]
  dzdz_theta = kappa * dz_theta.deriv('z')
  
  turbterm = dxdx_theta[0..-1,0..-1,2..99,0..-1] + dzdz_theta[2..1199,0..-1,0..-1,0..-1]
  
  turb_model = densBZ * turb

  turb_XYZmean = turb_model.mean('x').mean('y').mean('z')
  turb_XYmean = turb_model.mean('x').mean('y')
  turb_XYmeanZcut = turb_XYmean.cut('z'=>zcut)
p "turb"
p turb_XYZmean.cut('t'=>puttime).val.to_f

  #turb_XYZmean = turbterm.mean('x').mean('y').mean('z')
  #turb_cut = turb
  #turb_cut = turb[1..1199,0..-1,1..99,0..-1]
  #turb_XYZmean = turb_cut.mean('x').mean('y').mean('z')

# rhoQ

  rhoQ = densBZ * rad
  #rhoQ = densBZ[1..1199,0..-1,1..99] * rad[1..1199,0..-1,1..99,0..-1]
  
  rhoQ_XYZmean = rhoQ.mean('x').mean('y').mean('z')
  rhoQ_XYmean = rhoQ.mean('x').mean('y')
  rhoQ_XYmeanZcut = rhoQ_XYmean.cut('z'=>zcut)
p "rhoQ"
p rhoQ_XYZmean.cut('t'=>puttime).val.to_f

# flux term

  rho_theta_u = densBZ * ptemp * velx
  rho_theta_w = densBZ * ptemp * velz 
  
  dx_rho_theta_u = rho_theta_u.deriv('x')
  dz_rho_theta_w = rho_theta_w.deriv('z')
  
  fluxterm = dx_rho_theta_u[0..-1,0..-1,1..99,0..-1] + dz_rho_theta_w[1..1199,0..-1,0..-1,0..-1] 
  
  fluxterm_XYZmean = fluxterm.mean('x').mean('y').mean('z')
  fluxterm_XYmean = fluxterm.mean('x').mean('y')
  fluxterm_XYmeanZcut = fluxterm_XYmean.cut('z'=>zcut)

# rho adv
  adv_rho = densBZ * adv
  adv_XYZmean = adv_rho.mean('x').mean('y').mean('z')
  adv_XYmean = adv.mean('x').mean('y')
  adv_XYmeanZcut = adv_XYmean.cut('z'=>zcut)

  p adv_XYZmean.cut('t'=>cuttime).mean('t')


# +-
  diff = dt_rho_theta_XYZmean - adv_XYZmean[1..500] 
  diff_rad = rhoQ_XYZmean + turb_XYZmean
  diff_cut = dt_rho_theta_XYmeanZcut - ( adv_XYmeanZcut[1..500] + turb_XYmeanZcut[1..500] + rhoQ_XYmeanZcut[1..500])
#  all_cut = - adv_XYmeanZcut[1..400] + turb_XYmeanZcut[1..400] + rhoQ_XYmeanZcut[1..400]
  all_cut =  adv_XYmeanZcut + turb_XYmeanZcut + rhoQ_XYmeanZcut

DCL.gropn(ins)
#DCL.sgpset('lfull',true)
#var0 = var0[5..104,0..-1,5..104]
#VarAll = var0.cut(true,true,true,0.0) + var_base

#GGraph.line( conv_XYZmean.cut('t'=>cuttime), 
#             true, 'exchange'=>false ,
#             'index'=>22, 'type'=>1, 'label'=>'nonconv',
#             'title'=>'Mean Potemtial Temp.',
#             'max'=>comax, 'min'=>comin)    

#=begin
GGraph.line( dt_rho_theta_XYZmean.cut('t'=>cuttime), 
             true, 'exchange'=>false ,
             'index'=>12, 'type'=>1,# 'label'=>'rhotheta',
             'title'=>'PTemp Energy',
             'max'=>comax, 'min'=>comin)    
#=begin
GGraph.line( turb_XYZmean.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
             'index'=>32, 'type'=>1, 'label'=>'turb',
             'max'=>comax, 'min'=>comin)    

GGraph.line( rhoQ_XYZmean.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
             'index'=>42, 'type'=>1, 'label'=>'rhoQ',
             'max'=>comax, 'min'=>comin)    

#GGraph.line( fluxterm_XYZmean.cut('t'=>cuttime), 
#             false, 'exchange'=>false ,
#             'index'=>62, 'type'=>1, 'label'=>'flux',
#             'max'=>comax, 'min'=>comin)    
#=end
GGraph.line(  adv_XYZmean.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
             'index'=>22, 'type'=>1, #'label'=>'adv',
             'max'=>comax, 'min'=>comin)    
#GGraph.line( diff.cut('t'=>cuttime), 
#             false, 'exchange'=>false ,
#             'index'=>72, 'type'=>1, 'label'=>'diff',
#             'max'=>comax, 'min'=>comin)    
#GGraph.line( diff_rad.cut('t'=>cuttime), 
#             false, 'exchange'=>false ,
#             'index'=>52, 'type'=>1, 'label'=>'diff',
#             'max'=>comax, 'min'=>comin)    
#outdif = diff.cut('t'=>cuttime).mean('t')
#p outdif.val.to_f
#=end

=begin
GGraph.line( dt_rho_theta_XYmeanZcut.cut('t'=>cuttime), 
             true, 'exchange'=>false ,
             'index'=>12, 'type'=>1, 'label'=>'rhotheta',
             'title'=>"z=#{zcut}",
             'max'=>comax, 'min'=>comin)    
#=begin
GGraph.line( adv_XYmeanZcut.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
#              true, 'exchange'=>false ,
#             'title'=>"z=#{zcut}",
             'index'=>22, 'type'=>1, 'label'=>'adv',
             'max'=>comax, 'min'=>comin)    

#GGraph.line( diff_cut.cut('t'=>cuttime), 
#             false, 'exchange'=>false ,
#             'index'=>32, 'type'=>1, 'label'=>'dif',
#             'max'=>comax, 'min'=>comin)    
GGraph.line( turb_XYmeanZcut.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
             'index'=>42, 'type'=>1, 'label'=>'turb',
             'max'=>comax, 'min'=>comin)    
GGraph.line( rhoQ_XYmeanZcut.cut('t'=>cuttime), 
             false, 'exchange'=>false ,
             'index'=>52, 'type'=>1, 'label'=>'rhoQ',
             'max'=>comax, 'min'=>comin)    
#GGraph.line( all_cut.cut('t'=>cuttime), 
#             false, 'exchange'=>false ,
#             'index'=>62, 'type'=>1, 'label'=>'all',
#             'max'=>comax, 'min'=>comin)    
=end
DCL.grcls