#!/usr/bin/ruby

# Energy Balance Model

require "narray"
require './ebm_fig'
require './ebm'
require './ebm_param'


if ARGV.size == 0 then
  iws              = 1
  FigFN            = "dcl"
  sMode            = 1

  print "Select heat transfer type: \n"
  print "  1: Budyko-type\n"
  print "  2: Sellers-type\n"
  print "> "
  sIDHeatTransport = gets.to_i

  sSolarRadTOAMin  = 1100.0 / 4
  sSolarRadTOAMax  = 1900.0 / 4
  sAlbIceFree      = AlbedoIceFree
  sAlbIce          = AlbedoIce
  sTransferCoefC   = TransferCoefC
  sDiffCoefD       = DiffCoefD
else
  iws              = 2
  FigFN            = ARGV[0]
  sMode            = ARGV[1].to_i
  sIDHeatTransport = ARGV[2].to_i
  sSolarRadTOAMin  = ARGV[3].to_f
  sSolarRadTOAMax  = ARGV[4].to_f
  sAlbIceFree      = ARGV[5].to_f
  sAlbIce          = ARGV[6].to_f
  sTransferCoefC   = ARGV[7].to_f
  sDiffCoefD       = ARGV[8].to_f
end


sFlagDraw = true
#FlagDraw = false

jmax       = 128
sNDivision =  101

ebm = EBM.new( jmax, sTransferCoefC, sDiffCoefD, sAlbIce, sAlbIceFree )


a_CoolSolarRadTOA    = NArray.sfloat(sNDivision).fill(-1.0)
a_CoolLatIce         = NArray.sfloat(sNDivision).fill(-1.0)
a_CoolGMTemp         = NArray.sfloat(sNDivision).fill(-1.0)
a_WarmSolarRadTOA    = NArray.sfloat(sNDivision).fill(-1.0)
a_WarmLatIce         = NArray.sfloat(sNDivision).fill(-1.0)
a_WarmGMTemp         = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat1SolarRadTOA = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat1LatIce      = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat1GMTemp      = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat2SolarRadTOA = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat2LatIce      = NArray.sfloat(sNDivision).fill(-1.0)
a_FixLat2GMTemp      = NArray.sfloat(sNDivision).fill(-1.0)

for l in 0..(sNDivision-1)
  # solar radiation flux per unit area
#  sSolarRadTOA    = 334.4
  sSolarRadTOA    = ( sSolarRadTOAMax - sSolarRadTOAMin ) / ( sNDivision - 1 ) * l + sSolarRadTOAMin


  ebm.ebm( sIDHeatTransport, sSolarRadTOA, 200.0 )
  a_CoolSolarRadTOA[l] = sSolarRadTOA * 4
  a_CoolLatIce     [l] = ebm.ice_line
  a_CoolGMTemp     [l] = ebm.get_gmtemp
#  draw_temp( ebmc.get_tempice, ebmc.get_lat, ebmc.get_temp )

  ebm.ebm( sIDHeatTransport, sSolarRadTOA, 400.0 )
  a_WarmSolarRadTOA[l] = sSolarRadTOA * 4
  a_WarmLatIce     [l] = ebm.ice_line
  a_WarmGMTemp     [l] = ebm.get_gmtemp


#  y_Lat   = ebmc.get_lat
#  draw_temp( ebmc.get_tempice, y_Lat, y_TempC, y_TempW ) if sFlagDraw

#    print "SolarRad: ", sSolarRadTOA, ", Ice edge latitude : ", ebmc.ice_line, ", ", ebmw.ice_line, "\n"
end

for l in 0..(sNDivision-1)
  sLatIce = ( 90.0 - 0.0 ) / ( sNDivision - 1 ) * l + 0.0

  sSolarRadTOA = ebm.ebm_fixlat( sIDHeatTransport, sLatIce, 1 )
  a_FixLat1SolarRadTOA[l] = sSolarRadTOA * 4
  a_FixLat1LatIce     [l] = ebm.ice_line
  a_FixLat1GMTemp     [l] = ebm.get_gmtemp

  sSolarRadTOA = ebm.ebm_fixlat( sIDHeatTransport, sLatIce, 2 )
  a_FixLat2SolarRadTOA[l] = sSolarRadTOA * 4
  a_FixLat2LatIce     [l] = ebm.ice_line
  a_FixLat2GMTemp     [l] = ebm.get_gmtemp

#    print "SolarRad: ", sSolarRadTOA, ", Ice edge latitude : ", ebml.ice_line, "\n"
end

draw_preparation( iws, FigFN ) if sFlagDraw

case sMode
when 1 then
  case sIDHeatTransport
  when 1 then
    title1 = sprintf( "Budyko-type, C = %#.2f (W m-2 K-1)", sTransferCoefC )
    title2 = sprintf( "As = %#.2f, Ai = %#.2f", sAlbIceFree, sAlbIce )
  when 2 then
    title1 = sprintf( "Sellers-type, D = %#.2f (W m-2 K-1)", sDiffCoefD )
    title2 = sprintf( "As = %#.2f, Ai = %#.2f", sAlbIceFree, sAlbIce )
  end
  draw_ice_lat( title1, title2, sSolarRadTOAMin * 4, sSolarRadTOAMax * 4, a_FixLat1SolarRadTOA, a_FixLat1LatIce, a_FixLat2SolarRadTOA, a_FixLat2LatIce, a_CoolSolarRadTOA, a_CoolLatIce, a_WarmSolarRadTOA, a_WarmLatIce )
when 2 then
  draw_gmtemp( sSolarRadTOAMin * 4, sSolarRadTOAMax * 4, 200.0, 400.0, a_FixLat1SolarRadTOA, a_FixLat1GMTemp, a_FixLat2SolarRadTOA, a_FixLat2GMTemp, a_CoolSolarRadTOA, a_CoolGMTemp, a_WarmSolarRadTOA, a_WarmGMTemp )
else
  print "Unexpected mode\n"
  exit
end

draw_finish  if sFlagDraw