#!/usr/bin/env ruby
# ----------------------------------------------
# local load path

$local_path = '/home/yukiko/work/ape/yukiko/lib'
$: << $local_path
$fig_path = "/home/yukiko/work/ape/NetCDF/other_group/figs/tmp/"

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

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


def set_directry
  $groupid = $groupid_hash[$rezol]
  $ncfile_path = "/home/yukiko/work/ape/NetCDF/other_group/#{$rezol.downcase}/"
  $file_label = "#{$rezol}_#{$expID}"
end

END{

#    sstid = ["control5n","1keq","3keq","3kw1","peaked","control","flat","qobs"]
    sstid = ["control"]
    sstid.each{ |item|
      $expID = item
      $rezol = "ECMWF_07"
      set_directry
      ecmwf_mk_pf

  }

}

##------

def ecmwf_mk_pf


  ofile = NetCDF.create("#{$groupid}_PFp_#{$expID}.nc")

  ifile = ape_new("#{$ncfile_path}#{$groupid}_PF_#{$expID}.nc")
  ps  = ape_new("#{$ncfile_path}#{$groupid}_SH_#{$expID}.nc").
    gphys_open("sh_ps")
#  plev = NArray[1000, 925, 850, 700, 600, 500, 400, 300, 250, 200, 
#    150, 100, 70, 50, 30, 20, 10]
  plev = NArray[1000, 2000, 3000, 5000, 7000, 10000, 15000, 20000, 25000, 30000, 40000, 50000, 60000, 70000, 85000, 92500, 100000]
  
  ["pf_q","pf_q_cld","pf_q_conv","pf_t","pf_t_cld","pf_t_conv","pf_t_lw","pf_t_sw","pf_u","pf_u_conv","pf_v","pf_v_conv"].each{ |var|
#  ["pf_q"].each{ |var|
    
    gphys = ifile.gphys_open(var)
    gphys_out = mks2p(gphys,ps,plev)
    gphys_out = gphys_out 
    GPhys::NetCDF_IO.write(ofile, gphys_out )
    
  }

  # 大域属性
  ofile.put_att("Conventions", "CF-1.0")
  ofile.put_att("history", "Original data produced: #{Time.now} from PF netCDF by Yukiko YAMADA (AGU for APE)")
  ofile.put_att("institution", "#{$rezol}")
  ofile.close

end


def mks2p(gphys,ps,plev)

#  garr1 , gphys(im,jm,km)
#  garr2 , gphys_out(im,jm,plev.size)
#  gp1 , plev_org(im,jm,km)
#  gp2 , plev(plev.size)

  plev_org = sigma_ps(ps)

  grid_i = gphys.grid_copy.axis(0)
  grid_j = gphys.grid_copy.axis(1)
  im = gphys.grid_copy.axis(0).pos.val.size
  jm = gphys.grid_copy.axis(1).pos.val.size
  km = gphys.grid_copy.axis(2).pos.val.size
  grid_k =
    Axis.new().
    set_pos(VArray.new(plev).rename("plev").
              put_att("units","hPa"))

  wkgp1 = NArray.sfloat(km).fill(0)
  wkgarr1 = NArray.sfloat(km).fill(0)
  gphys_out = NArray.sfloat(im,jm,plev.size).fill(0)

  jm.times {|j|
    im.times {|i|

      if plev_org[i,j,0] <= plev_org[i,j,1] 
        km.times {|k|
          wkgp1[ k ]   = plev_org[ i, j, k ]
          wkgarr1[ k ] = gphys[ i, j, k ,0].data.val.to_na
        }
      else
        km.times {|k|
          wkgp1[ k ]   = plev_org[ i, j, km-(k+1) ]
          wkgarr1[ k ] = gphys[ i, j, km-(k+1),0 ].data.val.to_na
        }
      end 
        
      plev.size.times { |k|
        
        if ( wkgp1[ 0 ] > plev[ k ] ) then
          gphys_out[ i, j, k ] = -999.0
        else          
          l = km
          1.upto(km-1) {|ll|
#            p "#{wkgp1[ ll ]}, #{plev[ k ]}"
            if ( wkgp1[ ll ] >= plev[ k ] ) then
              l = ll
              break
            end 
          }
          
          if( l == km ) then
            gphys_out[ i, j, k] = -999.0
          else
            gphys_out[ i, j, k ] = 
              ( wkgarr1[ l ] - wkgarr1[ l-1 ] )/( log( wkgp1[ l ]/wkgp1[ l-1 ] ) ) * ( log( plev[ k ] / wkgp1[ l-1 ] ) ) + wkgarr1[ l-1 ]

#             p gphys_out[ i, j, k ]
          end 
       end 
        }

#      p "#{i}, #{j} #{gphys_out[i,j,true].to_a}"
#      p "#{i}, #{j}"
}
}

  grid = Grid.new(grid_i,grid_j,grid_k)
  gphys_out = GPhys.new(grid, VArray.new(gphys_out).rename(gphys.name).
              put_att("units",gphys.data.get_att("units")).
              put_att("ape_name",gphys.data.get_att("ape_name")) )

return gphys_out

end


def sigma_ps(ps)

   a = NArray[ 0.0,
    20.000000,
    38.425343,
    63.647804,
    95.636963,
   134.483307,
   180.584351,
   234.779053,
   298.495789,
   373.971924,
   464.618134,
   575.651001,
   713.218079,
   883.660522,
  1094.834717,
  1356.474609,
  1680.640259,
  2082.273926,
  2579.888672,
  3196.421631,
  3960.291504,
  4906.708496,
  6018.019531,
  7306.631348,
  8765.053711,
 10376.126953,
 12077.446289,
 13775.325195,
 15379.805664,
 16819.474609,
 18045.183594,
 19027.695313,
 19755.109375,
 20222.205078,
 20429.863281,
 20384.480469,
 20097.402344,
 19584.330078,
 18864.750000,
 17961.357422,
 16899.468750,
 15706.447266,
 14411.124023,
 13043.218750,
 11632.758789,
 10209.500977,
  8802.356445,
  7438.803223,
  6144.314941,
  4941.778320,
  3850.913330,
  2887.696533,
  2063.779785,
  1385.912598,
   855.361755,
   467.333588,
   210.393890,
    65.889244,
     7.367743,
     0.000000,
     0.000000] 

  b = NArray[ 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000000000,
 0.0000758235,
 0.0004613950,
 0.0018151561,
 0.0050811190,
 0.0111429105,
 0.0206778757,
 0.0341211632,
 0.0516904071,
 0.0735338330,
 0.0996746942,
 0.1300225109,
 0.1643843204,
 0.2024759352,
 0.2439331412,
 0.2883229554,
 0.3351548910,
 0.3838921487,
 0.4339629412,
 0.4847715795,
 0.5357099175,
 0.5861684084,
 0.6355474591,
 0.6832686067,
 0.7287858129,
 0.7715966105,
 0.8112534285,
 0.8473749161,
 0.8796569109,
 0.9078838825,
 0.9319403172,
    0.9518215060,
    0.9676452279,
    0.9796627164,
    0.9882701039,
    0.9940194488,
    0.9976301193,
    1.0000000000 ]

  data = NArray.sfloat(ps[true,0,0].data.val.size ,ps[0,true,0].data.val.size, 47).fill(0)

  47.times{ |num|
    numnum = num + 14
  data[true,true,num] = a[numnum] + b[numnum]*ps[true,true,0].data.val.to_na
  }

  return data

end