#!/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{

$rezol = ["T39L48_eml",
  "T39L24_eml",
  "T39L96_eml",
  "T79L48_eml" ]

# $rezol = ["T79L48_non"]
#$rezol = ["T39L48_eml"]
#$resol = "T39L48_eml"
#$resol = "T159L48_non"
#$resol = "T319L48_non"
#$sstid = "control"

 sstid_all = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N","flat3keq","Qobs3keq","H1998con","H1998pa"]
 sstid = ["control","flat","peaked","1keq","3keq","3kw1","Qobs","control-5N"]
sstid_tmp = ["control","flat","Qobs","3keq","flat3keq","Qobs3keq",
  "H1998con","H1998pa"]

    $sstid = "control"
    $resol = "CSIRO"
    set_directry
    mknc_stspctrum

    $resol = "MGO"
    set_directry
    mknc_stspctrum

    $sstid = "controln48"
    $resol = "UKMO"
    set_directry
    mknc_stspctrum


}

def set_directry
  $groupid = $groupid_hash[$resol]
  $ncfile_path = "/work11/ape/NetCDF/other_group/#{$resol.downcase}/"
  $file_label = "#{$resol}_#{$expID}"
end

def mknc_stspctrum

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

  $t= ape_new "#{$ncfile_path}#{$groupid}_TR_#{$sstid}.nc"    
#  $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"

  $t.netcdf_open.var_names.each { |name| 
    if name =~ /lw_toa|mslp|500|250|ps|tppn/
    
      puts name
      gphys = $t.go(name)
#      gphys = __tr_T159_marge(name)
#      gphys = __tr_T319_marge(name)
      
#      gphys_sym, gphys_asym, gphys_sym_wk, gphys_asym_wk, gphys_bg = spct_wk1999_ichimai(gphys)
      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)
      
    end
  }      

  # 大域属性
  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 by Yukiko YAMADA (AGU for APE)")
  file.put_att("institution", "#{$resol}")
  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 spct_wk1999(gphys)
  
  gphys1 = gphys.cut(true,-15..0,true)
  gphys2 = gphys.cut(true,0..15,true)

#  gphys1 = gphys.cut(true,-5..0,true)
#  gphys2 = gphys.cut(true,0..5,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_wk1999_ichimai(gphys)
  
  gphys = gphys.cut(true,0,true)[true,-401..-1]
  gphys = gphys.stspct_fft("spct")

  gphys_sym  = gphys.
    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_bg  = spct_121mean(gphys_sym)

  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_bg = gphys_bg.rename("#{gphys.name}_bg_spct").
    set_att("ape_name", gphys.data.get_att("ape_name")).
    set_att("units", "1" )

  gphys_asym = 0.0
  gphys_asym_wk = 0.0

  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|
#  80.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