#!/bin/csh
#  D. Sandwell 1/18/03
#
#  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 < 8) then
  echo " Usage: update_grid topo.xyz topo_8.2.img lon0 lonf lat0 latf dlon dlat"
  echo "        topo.xyz   -  file of lon, lat, depth (neg meters)"
  echo "        topo_8.2.img - global 2-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 "  "
  echo " Note the input data may have negative longitude but the"
  echo " output grid will be longitude east.  Also this code will not"
  echo " span Greenwich."
  echo " "
  echo "Example: update_grid german_iceland.xyz topo_8.2.img 340 343 65 68 .01 .01"

  exit 1
 endif
#
#  interpolate the global grid
#
interp_ship_80 << END
$2
2
1
$1
1 2 3
temp.xyzz
END
#
#  dump the global grid data for the area
#
img2xyz_80 << END
$2
2
1
$5 $6 $3 $4
0
global.xyz
END
#
#  
#   make sure the longitudes are positive and then do a block median of the
#   difference.  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($1 < 0) $1 = $1 + 360; print ($1, $2, $4-$3)}' temp.xyzz > temp.xyd
awk '{print ($1, $2, "0.0")}' global.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 temp.xyzz temp.xyd block.xyd global.xyz global.grd diff.grd