#
# ximg01.rb
#
# $Id: ximg01.rb,v 1.2 2000/11/27 01:49:23 keiko Exp $
#

#-----------------------------------------------------------------------
#     TOMS OZONE DISTRIBUTION (MONTHLY MEAN CLIMATOLOGY)
#-----------------------------------------------------------------------

require "narray"
require "numru/dcl"

include NumRu
include Math


KXBK = 12
IXBK = KXBK*2+1
JXBK = 37
IXBG = 73
JXBG = 37
IXGX = 361
JXGX = 181
NCL = 100
NCL1 = NCL+1

bkdata = NArray.sfloat(IXBK*JXBK)
bgdata = NArray.sfloat(IXBG, JXBG)
gxdata = NArray.sfloat(IXGX, JXGX)
tlevn = NArray.sfloat(NCL1)

def fftezb(fc, nx)
  nn = 0
  maxwn = fc.shape[0]

  nmax = 512
  if (nx > nmax)
    cmsg = format("WORKING AREA IS NOT ENOUGH / NX SHOULD BE .LE. %4d.", nmax)
    DCL::msgdmp('E', 'FFTEZB', cmsg)
  end
  if (nx < maxwn*2)
    cmsg = "MAXWN CANNOT BE CALCULATED / NX SOULD BE .GE. MAXWN*2.0"
    DCL::msgdmp('E', 'FFTEZB', cmsg)
  end

  if (nx != nn)
    wsave = DCL::rffti(nx)
    nn=nx
  end
  y = Array.new(nx, 0)
  y[0] = fc[0]
  1.step(maxwn-1, 2) do |i|
    y[i] = fc[i] * 0.5
    y[i+1] = fc[i+1] * -0.5
  end
  y = DCL::rfftb(y, wsave)
end

def bktobg(bkdata, bgdata)
  rmiss = DCL::glrget('RMISS')
#     CALL DTPGET('MAXWN',KX)
  for j in 0..JXBG-1
    if (bkdata[0,j] == rmiss)
      bgdata[true,j] = rmiss
    else
      wk = fftezb(bkdata[true,j], IXBG-1)
      bgdata[0..-2, j] = wk
      bgdata[-1,j]=bgdata[0,j]
    end
  end
end

#-- data ---
ipatn = []
for i in 0..NCL-1
  ipatn[i] = ([11, i+1].max) *1000 + 999
end
tlevn.indgen(203.0, 3.0)

cdsn = 'tomsclm.dat'
file = open(cdsn, "r")

DCL::swlstx('LWAIT', false)
DCL::swlstx('LWAIT0', false)
DCL::swlstx('LDUMP', false)
DCL::swistx('IPOSX', 200)
DCL::swistx('IPOSY', 200)
DCL::swistx('IWIDTH', 800)
DCL::swistx('IHEIGHT', 400)

#-- graph ---
iws = (ARGV[0] || (puts ' WORKSTATION ID (I)  ? ;'; DCL::sgpwsn; gets)).to_i
DCL::gropn iws
DCL::gllset('LMISS', true)
rmiss = DCL::glrget('RMISS')
DCL::sglset('LFULL', true)
DCL::sglset('LSOFTF', false)
DCL::umlset('LGLOBE', true)
DCL::slrat(1.0, 0.5)

for m in 1..12
  mon = file.gets.split(" ")[-1].to_i
  bkdata.reshape(IXBK*JXBK)
  i = 0
  (JXBK*3).times do |j|
    file.gets.split(" ").each do |data|
      bkdata[i] = data.to_f
      i += 1
    end
  end
  bktobg(bkdata.reshape(IXBK, JXBK), bgdata)
  gxdata.fill!(rmiss)
  for j in 0..JXBG-1
    for i in 0..IXBG-1
      gxdata[i*5, j*5] = bgdata[i,j]
    end
    wk = DCL::vrintr(gxdata[true, j*5], IXGX, 1)
    gxdata[true, j*5] = wk
  end
  for i in 0..IXGX-1
    wk = DCL::vrintr(gxdata[i, true], JXGX, 1)
    gxdata[i, true] = wk
  end

  DCL::grfrm
  DCL::grswnd(0.0, 360.0, -90.0, 90.0)
  DCL::grsmpl(180.0, 90.0, 0.0)
  DCL::grstrn(DCL::isgtrc('HMR'))
  DCL::umpfit
  DCL::grstrf

  DCL::slpwwr(1)
  DCL::ueitlv
  DCL::uestln(tlevn, ipatn)
  DCL::uetonf(gxdata)

  DCL::umpgrd
  DCL::umplim
  DCL::umpmap('coast_world')
end
file.close

DCL::grcls