** use 3-day running mean for IR image
function plot(args)
time1=subwrd(args, 1)
*time1=12Z27nov2000
say time1

'reinit'
'open /export-4/sgi109/pacs/mcidas/data/global_ir/global_ir_geos_merged_half-deg/global_0.5deg_0.5hourly_tmp-lnx.ctl'
'open daily.clim.200chi.71-00.harm4.ctl'
'open /export-5/cacsrv1/gdas/6hr/pgb.ctl'
'enable print /disk2/irstat/web_plots/meta/am_ir_monthly.meta'
'set gxout shaded'
'set grid off'
*'set mpdset mres'
*'run bw.gs'
*'run /export/lnx196/yzhou/CWlink/IR/rgbset.gs'
'set grads off'
'set map 0 1 1'
'set lat -60 60'
'set lon -180 180'
'set time 'time1
'q dims'
line5=sublin(result,5)
tmp2=subwrd(line5,9)
'set t 'tmp2-72
'draw string 5.15 1.47 Degrees K'
'run /disk2/irstat/web_plots/web.specs.monthly'
say 'd ave(tmpm,t='tmp2-143',t='tmp2',1,-b)/10'
'd ave(tmpm,t='tmp2-143',t='tmp2',1,-b)/10'
*'d ave(tmpm,t-1,t+0)/10'
'run /disk2/irstat/web_plots/cbar1.gs'
*'run /export/lnx196/yzhou/CWlink/IR/date.gs'
'draw string 5.5 6.55 'mydate()
say mydate()

** compute velocity potential anom from current value and clim
'set dfile 3'
'set time 'time1
'set lat -60 60'
'set lon -180 180'
'set lev 200'
say 'define vpot1=ave(vpot.3,t-11,t+0,1,-b)'
'define vpot1=ave(vpot.3,t-11,t+0,1,-b)'

** year in control file for climo should be the current year
'set dfile 2'
'set t last'
'q dims'
line5=sublin(result,5)
ndays=subwrd(line5,9)

'set time 'time1
'q dims'
line5=sublin(result,5)
tmp2=subwrd(line5,9)

if (tmp2 > 1 ) 
tmp1=tmp2-1
else
tmp1=ndays
endif

if (tmp2=ndays)
tmp3=1
else
tmp3=tmp2+1
endif

'set lat -60 60'
'set lon -180 180'
'set lev 200'
say 'define vpot0=(vpot.2(t='tmp1')+vpot.2(t='tmp2')+vpot.2(t='tmp3'))/3.0'
'define vpot0=(vpot.2(t='tmp1')+vpot.2(t='tmp2')+vpot.2(t='tmp3'))/3.0'

say 'define anom=vpot1-vpot0'
'run rgbset.gs'
'define anom=vpot1-vpot0'
'set gxout contour'
'set clab off'
'set ccolor 1'
'set cthick 8'
'set grads off'
'set clevs  -12e+6 -10e+6 -8e+6 -6e+6 -4e+6 -2e+6 2e+6  4e+6 6e+6 8e+6 10e+6 12e+6' 
'set ccols 36 36 36 36 36 36 77 77 77 77 77 77' 
*'set cstyle 2 2 2 2 2 2 1 1 1 1 1 1'
'd anom'

'print'
'disable print'
'quit'

function mydate
  'query time'
  sres = subwrd(result,3)
  i = 1
  while (substr(sres,i,1)!='Z')
    i = i + 1
  endwhile
  hour = substr(sres,1,i)
  isav = i
  i = i + 1
  while (substr(sres,i,1)>='0' & substr(sres,i,1)<='9')
    i = i + 1
  endwhile
  day = substr(sres,isav+1,i-isav-1)
  month = substr(sres,i,3)
  year = substr(sres,i+3,4)
  return (day' 'month' 'year)