require 'getopts'
require 'numru/gphys'

module NumRu

  class NetCDFVar

    def get_with_miss_and_scaling2(*args) # make mask before scaling
      __interpret_missing_params if !defined?(@missval)
      packed_data = simple_get(*args)
      scaled_data = scaled_get(*args)
      sf = att('scale_factor')
      ao = att('add_offset')
      if @vmin || @vmax
	if sf && ao
	  csf = sf.get
	  cao = ao.get
	  vmin = (@vmin-cao)/csf if @vmin
	  vmax = (@vmax-cao)/csf if @vmax
	elsif
	  vmin = @vmin; vmax = @vmax
	end
	if vmin
	  mask = (packed_data >= vmin) 
	mask = mask.and(packed_data <= vmax) if vmax
	else
	  mask = (packed_data <= vmax)
	end
	data = NArrayMiss.to_nam(scaled_data, mask)
      elsif @missval	# only missing_value is present.
	eps = 1e-6
	missval = @missval[0].to_f
	vmin = missval - eps
	vmax = missval + eps
	mask = (packed_data <= vmin) 
	mask = mask.or(packed_data >= vmax)
	data = NArrayMiss.to_nam(scaled_data, mask)
      else
	data = scaled_data
      end
      data
    end

    def put_with_miss_after_scaling(data, *args)
      if data.is_a?( NArrayMiss )
	__interpret_missing_params if !defined?(@missval)
	if @missval
	  scaled_put_without_missval(data, *args)
	else
	  scaled_put(data.to_na, *args)
	end
      else
	scaled_put(data, *args)
      end
    end

    def scaled_put_without_missval(var,hash=nil)
      sf = att('scale_factor')
      ao = att('add_offset')
      if ( sf == nil && ao == nil ) then
	# no scaling --> just call put
	simple_put(var,hash)
      else
	if (sf != nil) 
	  csf = sf.get
	  if csf.is_a?(NArray) then  # --> should be a numeric
	    csf = csf[0]
	  elsif csf.is_a?(String)
	    raise TypeError, "scale_factor is not a numeric"
	  end
	  if(csf == 0) then; raise NetcdfError, "zero scale_factor"; end
	else
	  csf = 1.0      # assume 1 if not defined
	end
	if (ao != nil) 
	  cao = ao.get
	if cao.is_a?(NArray) then  # --> should be a numeric
	  cao = cao[0]
	elsif csf.is_a?(String)
	  raise NetcdfError, "add_offset is not a numeric"
	end
	else
	  cao = 0.0      # assume 0 if not defined
	end
	var = var.to_na( @missval[0]*csf + cao)
	simple_put( (var-cao)/csf, hash )
      end
    end  
    
  end
end

################################################################################
include NumRu
include Misc::EMath

unless getopts("v:", "shum:", "mask:", "output:epflx.nc")
  print "#{$0}:illegal options.\n"
  exit 1
end

L   = UNumeric.new(2.5008e6, 'J.kg-1')

z_axs = 'level'
z_range = [300, 1000]

nc_vwnd, var_vwnd =   ($OPT_v).split(/\s*@\s*/)
nc_shum, var_shum =   ($OPT_shum).split(/\s*@\s*/)
nc_mask, var_mask =   ($OPT_mask).split(/\s*@\s*/)

gp_v =     GPhys::IO.open( nc_vwnd,  var_vwnd ).cut(z_axs=>z_range[0]..z_range[-1])
gp_shum =  GPhys::IO.open( nc_shum,  var_shum )
gp_mask =  GPhys::IO.open( nc_mask,  var_mask ).cut(z_axs=>z_range[0]..z_range[-1])

ofile = NetCDF.create($OPT_output)

nt = gp_v.shape[-1]
i = 0
GPhys::IO.each_along_dims_write([gp_v, gp_shum, gp_mask], ofile, -1){
  |v, shum, mask|
  i += 1
  print "processing #{i} / #{nt} ..\n" if (i % (nt/20+1))==1
  
  v    = v + mask
  shum = shum + mask

  [v, shum].each{|gp|                            #  This part will not
    gp.data.att_names.each{|nm|                  #  be needed in future.
      gp.data.del_att(nm) if /^valid_/ =~ nm     #  (Even now, it is not
      gp.data.del_att(nm) if /^title/  =~ nm     #  needed if the valid	 
    }                                            #  range is wide enough)
  }                                              

  v_bar     = v.mean('lon')
  v_dash    = v - v_bar
  shum_bar  = shum.mean('lon')
  shum_dash = shum - shum_bar
  
  lhefl = ( v * shum ).mean('lon') * L
  lhefl_bar  = ( v_bar  * shum_bar  ) * L
  lhefl_dash = ( v_dash * shum_dash ).mean('lon') * L 
  
  lhefl.data.rename!('lhefl')
  lhefl.data.set_att('long_name', 'latent heat energy transport')
  lhefl.data.set_att('units', 'W')
  lhefl_bar.data.rename!('lhefl_bar')
  lhefl_bar.data.set_att('long_name', 'latent heat energy transport (bar)')
  lhefl_bar.data.set_att('units', 'W')
  lhefl_dash.data.rename!('lhefl_dash')
  lhefl_dash.data.set_att('long_name', 'latent heat energy transport (dash)')
  lhefl_dash.data.set_att('units', 'W')
  
  
  [ lhefl, lhefl_bar, lhefl_dash ]
}

  ofile.close