function composite()

* Fresh start
* -----------
  'reinit'

* Add contribution of each site
* -----------------------------
  accumsite(abr,1)
  sites1()
  sites2()
*  sites3()

*  Must have at least 1 file open
*  -------------------------------
   opensite(abr)
   'set lon -135 -64'
   'set lat 20 55'

* Finally, normalize to avoid double counts
* -----------------------------------------
  'define vc = if(nn,>,0,vv/nn,-u)'

* Save the composite to a netcdf file
* -----------------------------------
  'define cvalue = if(nn,>,0,vv/nn,vv)'
  'set sdfwrite radar_composite.nc'
  'sdfwrite cvalue'

* Simple plot
* -----------
  'set mpdset hires'
  'set gxout grfill'
  'd vc'
  'draw title Current NEXRAD Radar'
  'cbarn'
* 'gxyat -x 20000 -y 10200 radar.png' 
  'gxyat -x 1600 -y 1200 radar.png' 



function sites1()
  accumsite(abx)
  accumsite(akq)
  accumsite(ama)
  accumsite(amx)

function sites2()
accumsite(apx)
accumsite(arx)
accumsite(atx)
accumsite(bbx)
accumsite(bgm)
accumsite(bhx)
accumsite(bis)
accumsite(blx)
accumsite(bmx)
accumsite(box)
accumsite(bro)
accumsite(buf)
accumsite(byx)
accumsite(cae)
accumsite(cbw)
accumsite(cbx)
accumsite(ccx)
accumsite(cle)
accumsite(clx)
accumsite(crp)
accumsite(cxx)
accumsite(cys)
accumsite(dax)
accumsite(ddc)
accumsite(dfx)
accumsite(dgx)
accumsite(dix)
accumsite(dlh)
accumsite(dmx)
accumsite(dox)
accumsite(dtx)
accumsite(dvn)
accumsite(dyx)
accumsite(eax)
accumsite(emx)
accumsite(enx)
accumsite(eox)
accumsite(epz)
accumsite(esx)
accumsite(evx)
accumsite(ewx)
accumsite(eyx)
accumsite(fcx)
accumsite(fdr)
accumsite(fdx)
accumsite(ffc)
accumsite(fsd)
accumsite(fsx)
accumsite(ftg)
accumsite(fws)
accumsite(ggw)
accumsite(gjx)
accumsite(gld)
accumsite(grb)
accumsite(grk)
accumsite(grr)
accumsite(gsp)
accumsite(gwx)
accumsite(gyx)
accumsite(hdx)
accumsite(hgx)
accumsite(hnx)
accumsite(hpx)
accumsite(htx)
accumsite(ict)
accumsite(icx)
accumsite(iln)
accumsite(ilx)
accumsite(ind)
accumsite(inx)
accumsite(iwa)
accumsite(iwx)
accumsite(jax)
accumsite(jgx)
accumsite(jkl)
accumsite(lbb)
accumsite(lch)
accumsite(lix)
accumsite(lnx)
accumsite(lot)
accumsite(lrx)
accumsite(lsx)
accumsite(ltx)
accumsite(lvx)
accumsite(lwx)
accumsite(lzk)
accumsite(maf)
accumsite(max)
accumsite(mbx)
accumsite(mhx)
accumsite(mkx)
accumsite(mlb)
accumsite(mob)
accumsite(mpx)
accumsite(mqt)
accumsite(mrx)
accumsite(msx)
accumsite(mtx)
accumsite(mux)
accumsite(mvx)
accumsite(mxx)
accumsite(nkx)
accumsite(nqa)
accumsite(oax)
accumsite(ohx)
accumsite(okx)
accumsite(otx)
accumsite(pah)
return

function sites3()
accumsite(pbz)
accumsite(pdt)
accumsite(poe)
accumsite(pux)
accumsite(rax)
accumsite(rgx)
accumsite(riw)
accumsite(rlx)
accumsite(rtx)
accumsite(sfx)
accumsite(sgf)
accumsite(shv)
accumsite(sjt)
accumsite(sox)
accumsite(srx)
accumsite(tbw)
accumsite(tfx)
accumsite(tlh)
accumsite(tlx)
accumsite(twx)
accumsite(tyx)
accumsite(udx)
accumsite(uex)
accumsite(vax)
accumsite(vbx)
accumsite(vnx)
accumsite(vtx)
accumsite(vwx)
accumsite(yux)
return

#-----------------------------------------

function accumsite(site,init)

* Open a NetCDF file of radar data
* --------------------------------
  opensite(site)

* Set domain of interest
* ----------------------
  'set lon -135 -64'
  'set lat 20 55'

* value on reference grid
* -----------------------
***'define v = re(value,0.013599)'
'define v = re(value,0.05)'

* has 1 on gridpoints with a value, zero elsewhere
* ------------------------------------------------
'define n = const(0*v+1,0,-u)'

* Now accumulate
* --------------
if (init=1) 
    'define vv = const(v,0,-u)'
    'define nn = n'
else
    'define vv = vv + const(v,0,-u)'
    'define nn = nn + n'
endif

* close the file
* --------------
  'close 1'

return

function opensite(site)
   say 'working on out/sn.k' site '.nc'
   'sdfopen out/sn.k' site '.nc'
return