#!/usr/bin/env ruby

# ----------------------------------------------
# local load path

$local_path = '/work11/ape/yukiko/lib' 
$: << $local_path

# ----------------------------------------------
# 必要なライブラリ, モジュールの読み込み

load "#{$local_path}/ape-view.rb"
load "#{$local_path}/lib-ape-view.rb"
load "#{$local_path}/lib-ape-base.rb"
load "#{$local_path}/lib-ape-stspctrum.rb"

# --------------------------------------------------------------------


END{

  $resol = "T159L48_eml"
  $sstid = "control"
  mknc_stspctrum


}

def mknc_stspctrum

  $ncfile_path = "/work11/ape/NetCDF/#{$resol}/"
  $groupid = "AGUforAPE-03a"
  $expID = $sstid

    
#  $t= ape_new "/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}.nc"    
#  $t= ape_new "/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half1.nc"    
#  $t= ape_new "/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_0031.nc"
  file = NetCDF.create("#{$resol}_#{$sstid}_wk_40smooth.nc")
  
  puts "#{$resol}_#{$sstid}_wk_40smooth.nc create start"

    
      puts "tr_tppn"
#      gphys = $t.go(name)
#      gphys = __tr_T159_marge(name)
#      gphys = __tr_T319_marge(name)

  gphys = __tr_T159_eml_marge("tr_tppn")

  gphys_sym, gphys_asym, gphys_sym_wk, gphys_asym_wk, gphys_bg = spct_wk1999(gphys)
  dim1 = gphys_sym.coord(0).shape.to_s.to_i
  dim2 = gphys_sym.coord(1).shape.to_s.to_i
      
  # netCDF 初期化
  GPhys::NetCDF_IO.write(file, gphys_sym_wk )
  GPhys::NetCDF_IO.write(file, gphys_asym_wk )
  GPhys::NetCDF_IO.write(file, gphys_bg)
  GPhys::NetCDF_IO.write(file, gphys_sym)
  GPhys::NetCDF_IO.write(file, gphys_asym)
      
  # 大域属性
  file.put_att("Conventions", "CF-1.0")
  file.put_att("title","Aqua Planet: #{$resol}_#{$sstid} Experiment, Wheeler and Kiladis, 1999 Plot")
  file.put_att("history", "Original data produced: #{Time.now} from TR netCDF")
  file.put_att("institution", "AGU for APE")
  file.put_att("source","afes_115_ape_#{$cumulus}.LM")
  file.close
  
end

# --------------------------------------------------------------------
  
def __tr_T159_marge(name)
  
  file1 = Ape.new("/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half1.nc")
  file2 = Ape.new("/work11/ape/NetCDF/#{$resol}/AGUforAPE-03a_TR_#{$sstid}_half2.nc")

  # data
  data = NArray.sfloat(480,40,1440).indgen * 0.0
  
  gphys = file1.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,0..719] = gphys.data.val
  
  gphys = file2.gphys_open(name)
  gphys = gphys.cut(true,-15..15,true)
  data[true,true,720..-1] = gphys.data.val
  
  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))


  # 軸
#  x_size = NArray.sfloat(480).indgen*360.0/480
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end

# --------------------------------------------------------------------
  
def __tr_T319_marge(name)
   

  file1 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0031.nc")
  file2 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0032.nc")
  file3 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0033.nc")
  file4 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0034.nc")
  file5 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0035.nc")    
  file6 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0036.nc")
  file7 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0037.nc")
  file8 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0038.nc")
  file9 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0039.nc")
  file10 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0040.nc")
  file11 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0041.nc")
  file12 = Ape.new("#{$ncfile_path}#{$groupid}_TR_#{$expID}_0042.nc")
  
    
  # data
  data = NArray.sfloat(960,80,1440).indgen * 0.0
  
  
  gphys = file1.gphys_open(name).cut(true,-15..15,true)
  data[true,true,0..119] = gphys.data.val
  gphys = file2.gphys_open(name).cut(true,-15..15,true)
  data[true,true,120..239] = gphys.data.val
  gphys = file3.gphys_open(name).cut(true,-15..15,true)
    data[true,true,240..359] = gphys.data.val
  gphys = file4.gphys_open(name).cut(true,-15..15,true)
  data[true,true,360..479] = gphys.data.val
  gphys = file5.gphys_open(name).cut(true,-15..15,true)
  data[true,true,480..599] = gphys.data.val
  gphys = file6.gphys_open(name).cut(true,-15..15,true)
  data[true,true,600..719] = gphys.data.val
  gphys = file7.gphys_open(name).cut(true,-15..15,true)
  data[true,true,720..839] = gphys.data.val
  gphys = file8.gphys_open(name).cut(true,-15..15,true)
  data[true,true,840..959] = gphys.data.val
  gphys = file9.gphys_open(name).cut(true,-15..15,true)
  data[true,true,960..1079] = gphys.data.val
  gphys = file10.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1080..1199] = gphys.data.val
  gphys = file11.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1200..1319] = gphys.data.val
  gphys = file12.gphys_open(name).cut(true,-15..15,true)
  data[true,true,1320..1439] = gphys.data.val
  
  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))
  
  # 軸
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end


# --------------------------------------------------------------------
  
def __tr_T159_eml_marge(name)

  $grfile_path = "/work11/ape/GrADS/T159L48_eml_expID01"

  # data
  data = NArray.sfloat(480,40,1440).indgen * 0.0

  num = "0011"
  gphys = Ape.new("#{$grfile_path}/#{num}/PRCPC.ctl"). 
    gphys_open("PRCPC").cut(true,-15..15,1,true)
  gphys_add = Ape.new("#{$grfile_path}/#{num}/PRCPL.ctl").
    gphys_open("PRCPL").cut(true,-15..15,1,true)
  gphys = gphys + gphys_add
  data[true,true,0..359] = gphys.data.val

  num = "0012"
  gphys = Ape.new("#{$grfile_path}/#{num}/PRCPC.ctl"). 
    gphys_open("PRCPC").cut(true,-15..15,1,true)
  gphys_add = Ape.new("#{$grfile_path}/#{num}/PRCPL.ctl").
    gphys_open("PRCPL").cut(true,-15..15,1,true)
  gphys = gphys + gphys_add
  data[true,true,360..719] = gphys.data.val

  num = "0013"
  gphys = Ape.new("#{$grfile_path}/#{num}/PRCPC.ctl"). 
    gphys_open("PRCPC").cut(true,-15..15,1,true)
  gphys_add = Ape.new("#{$grfile_path}/#{num}/PRCPL.ctl").
    gphys_open("PRCPL").cut(true,-15..15,1,true)
  gphys = gphys + gphys_add
  data[true,true,720..1079] = gphys.data.val

  num = "0014"
  gphys = Ape.new("#{$grfile_path}/#{num}/PRCPC.ctl"). 
    gphys_open("PRCPC").cut(true,-15..15,1,true)
  gphys_add = Ape.new("#{$grfile_path}/#{num}/PRCPL.ctl").
    gphys_open("PRCPL").cut(true,-15..15,1,true)
  gphys = gphys + gphys_add
  data[true,true,1080..-1] = gphys.data.val
    
  data = VArray.new(data).rename("#{name}").
    put_att("units",gphys.data.get_att("units")).
    put_att("ape_name",gphys.data.get_att("ape_name"))

  # 軸
  t_size = NArray.sfloat(1440).indgen*1.0/4.0 + 0.25
  y_size = gphys.grid_copy.axis(1).pos.val
  x_size = gphys.grid_copy.axis(0).pos.val
  
  x = VArray.new(x_size).rename("lon").put_att("units","degrees_east")
  y = VArray.new(y_size).rename("lat").put_att("units","degrees_north")
  t = VArray.new(t_size).rename("time").put_att("units","days since 0000-01-01")
  x = Axis.new().set_pos(x)
  y = Axis.new().set_pos(y)
  t = Axis.new().set_pos(t)
  grid = Grid.new(x,y,t)

  data = GPhys.new( grid, data)
  
  return data
end

# --------------------------------------------------------------------

def spct_wk1999(gphys)
  
  gphys1 = gphys.cut(true,-15..0,true)
  gphys2 = gphys.cut(true,0..15,true)

#  gphys1 = gphys.cut(true,-2..0,true)
#  gphys2 = gphys.cut(true,0..2,true)

  dim = gphys1.coord(1).shape.to_s.to_i  
  work1 = gphys1[true,0,0..359].stspct_fft("spct") * 0.0
  work2 = gphys2[true,0,0..359].stspct_fft("spct") * 0.0

  dim.times { |dim_num|
    10.times { |num|
      work1 = gphys1[true,dim_num,(num*120)..359+(num*120)].
	stspct_fft("spct") + work1
      work2 = gphys2[true,dim_num,(num*120)..359+(num*120)].
	stspct_fft("spct") + work2
    }
  }

  gphys_asym = (( work1 - work2 )/2 ).
    set_att("ape_name", "#{gphys.data.get_att("ape_name")} S-T power spectrum").
    set_att("units", "1" ). 
    rename("#{gphys.name}_sym_org_spct")
  gphys_sym  = (( work1 + work2 )/2 ).
    set_att("ape_name", "#{gphys.data.get_att("ape_name")} S-T power spectrum").
    set_att("units", "1" ). 
    rename("#{gphys.name}_asym_org_spct")
  gphys = ( gphys_sym + gphys_asym ).rename("#{gphys.name}")
  gphys_bg  = spct_121mean(gphys)

  gphys_sym_wk = spct_divide(gphys_sym,gphys_bg).
    set_att("ape_name", gphys.data.get_att("ape_name")).
    set_att("units","1"). 
    rename("#{gphys.name}_sym_spct")
  gphys_asym_wk = spct_divide(gphys_asym,gphys_bg).
    set_att("ape_name", gphys.data.get_att("ape_name") ).
    set_att("units","1"). 
    rename("#{gphys.name}_asym_spct")

  gphys_bg = gphys_bg.rename("#{gphys.name}_bg_spct").
    set_att("ape_name", gphys.data.get_att("ape_name")).
    set_att("units", "1" )

  return gphys_sym, gphys_asym, gphys_sym_wk, gphys_asym_wk, gphys_bg
end


def spct_121mean(gphys)
  
  dim1 = gphys.coord(0).shape.to_s.to_i
  dim2 = gphys.coord(1).shape.to_s.to_i
  
  gphys_bg = gphys.copy
  work     = gphys.copy

  40.times{ |num|
#  10.times{ |num|
    dim1.times { |dim1_num|
      if dim1_num == 0
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num,true] + gphys_bg[dim1_num+1,true] )/2 
      elsif dim1_num == (dim1-1)
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num-1,true] + gphys_bg[dim1_num,true] )/2 
      else
	work[dim1_num,true] = 
	  ( gphys_bg[dim1_num-1,true]+ gphys_bg[dim1_num,true] + gphys_bg[dim1_num,true] + gphys_bg[dim1_num+1,true] )/4 
      end
    }
    
    dim2.times { |dim2_num|
      if dim2_num == 0
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num] + work[true,dim2_num+1] )/2 
      elsif dim2_num == (dim2-1)
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num-1] + work[true,dim2_num] )/2 
      else
	gphys_bg[true, dim2_num] = 
	  ( work[true,dim2_num-1]+ work[true,dim2_num] + work[true,dim2_num] +  work[true,dim2_num+1] )/4 
      end
    }
  }    


  return gphys_bg

end

def spct_divide(gphys,gphys_bg)  

  dim1 = gphys.coord(0).shape.to_s.to_i
  dim2 = gphys.coord(1).shape.to_s.to_i

  dim1.times{ |dim1_num|
    dim2.times{ |dim2_num|
      unless gphys_bg[dim1_num,dim2_num].data.val == 0.0
	gphys[dim1_num,dim2_num]=
	  gphys[dim1_num,dim2_num].data.val/gphys_bg[dim1_num,dim2_num].data.val
      else 
	gphys[dim1_num,dim2_num] = 0.0 
      end
    }
  }

  return gphys
end




=begin
    $t.gropn(1)

    $t.plot(gphys_sym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    file_name = "#{$file_label}_#{gphys_sym.name}.gif"
    p $file_name

    # ダンプする場合
    GGraph.tone( gphys_sym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm` 

    $t.plot(gphys_asym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
    file_name = "#{$file_label}_#{gphys_asym.name}.gif"
    p $file_name
    
    # ダンプする場合
    GGraph.tone( gphys_asym[((dim1+1)/2-31)..((dim1+1)/2+29),0..80])
   `for file in *.xwd ; do xwdtopnm ${file} | pnmcut 2 2 904 654 > tmp.pnm; done`
    `ppmtogif tmp.pnm > #{$fig_path}#{file_name}`
    `rm dcl_* tmp.pnm` 

  }
  
  # おしまい
  $t.grcls
  `rm tmp.pnm dcl_*`
=end

=begin
  dim.times { |dim_num|
    5.times { |num|
      work1 = gphys1[true,dim_num,(num*240)..359+(num*240)].
	stspct_fft("spct") + work1
      work2 = gphys2[true,dim_num,(num*240)..359+(num*240)].
	stspct_fft("spct") + work2
    }
  }
=end