#!/bin/csh
#  D. Sandwell 08/08/07
#
#  Update the global topography grid using some new trusted data.
#  Note a similar program could easily be written to update the
#  gravity grid.
#
# check the number of arguments
#
 if ($#argv < 9) then
  echo " Usage: update_grid topo.xyz topo_10.1.img lon0 lonf lat0 latf dlon dlat edit"
  echo "        topo.xyz   -  file of lon, lat, depth (neg meters)"
  echo "        topo_10.1.img - global 1-min topography grid"
  echo "        lon0       -  left boundary of output grid (0-360)" 
  echo "        lonf       -  right boundary of output grid (0-360)"
  echo "        lat0       -  lower boundary of output grid"
  echo "        latf       -  upper boundary of output grid"
  echo "        dlon       -  longitude spacing of output grid (e.g., .001)"
  echo "        dlat       -  latitude spacing of output grid"
  echo "        edit       _  remove |differences| > edit"
  echo " "
  echo "Example: update_grid german_iceland.xyz topo_9.1.img 340 343 65 68 .01 .01"

  exit 1
 endif
echo $9
#
#  extract the global topo into a grd file as well as the
#  xyz constraint locations
#
#
img2grd $2 -T2 -S1 -V -R$3/$4/$5/$6 -m1 -D  -Gtrack.grd
grd2xyz track.grd -S > track.xyz
#
img2grd $2 -T1 -S1 -V -R$3/$4/$5/$6 -m1 -D -Gpredict.grd
grd2xyz predict.grd -S > global.xyz
#
#  interpolate the new data through the grid
#
grdtrack $1 -Gpredict.grd -S > predict.xyzz
#
#   We attached a bunch of zeros to the difference file at
#   the locations of the global grid.  This was done to make sure the 
#   surface program behaves well outside of the area constrained by real
#   data.  One could be more careful here to not include zeros in the area
#   of the high-density data but blockmedian may filter them anyway.
#
awk '{if(!(($4-$3) > '$9' || ($3-$4) > '$9')) print($1, $2, $4-$3)}' < predict.xyzz > temp.xyd
awk '{print ($1, $2, "0.0")}' track.xyz >> temp.xyd
blockmedian temp.xyd -R$3/$4/$5/$6 -I$7/$8 > block.xyd
#
#   now make the difference grid
#
surface block.xyd -R$3/$4/$5/$6 -I$7/$8 -T.75 -V -Gdiff.grd
#
#   make a matching grid from the global topo
#
surface global.xyz -R$3/$4/$5/$6 -I$7/$8 -T.75 -V -Gglobal.grd
#
#   add the two grids
#
grdmath  global.grd diff.grd ADD = final_topo.grd
#
#  clean up the mess
#
rm track.grd predict.grd global.xyz predict.xyz block.xyd global.grd temp.xyd
