#!/bin/bash
#
#   Example script to prepare the data and force grids for the itration flexure program
#   files in Mercator coordinates will be called _xy.grd
#   in mercator coordinates
#
#   first set the boundaries of the area
# 
    export BOUNDS="263/297/-40/-27.5"
    echo $BOUNDS
    gmtset ELLIPSOID = Sphere
#
#   get the gravity and topography grids using img2grd
#   and the age grid using grdcut
#
   ln -s /topex/ftp/pub/global_grav_1min/grav.img.23.1 .
   ln -s /topex/ftp/pub/global_topo_1min/topo_18.1.img .
   img2grd grav.img.23.1 -T1 -M -m1 -D -S.1 -R$BOUNDS -Ggrav_merc.grd -V
   img2grd topo_18.1.img -T2 -M -m1 -D -R$BOUNDS -Gtopo_merc.grd -V 
   img2grd topo_18.1.img -T1 -M -m1 -D -R$BOUNDS -Galltopo_merc.grd -V 

   ln -s ../../flex/age/age.cpt

#
# 	reproject the age grid to match the topo grid
#

#	adjusted bounds:
#	-R263/297/-40.0039553203/-27.4961156809
#

   grdproject age.grd -Jm111195 -N2040/906 -R$bounds -F -Sl -Gage_merc.grd

#
#  here is the output size
#  output will be 702 by 702 pixels.
#
#  next modify vload to match this size
#  edit vload.f and set
#  parameter(ni=718,nj=1200,nwork=32768,nj2=nj/2+1,ni2=ni/2+1)
#
#  kml2gmt sam_100k_line.kml | tail -n +4 > 100k_len.dat
   
	sed '$d'  cchile_100k.dat > one.dat
	tail -n +2 cchile_100k.dat > two.dat
	paste one.dat two.dat | awk '{print $1,$3,$2,$4}' > load_segs.dat
    awk '{print $1,$2,$3,$4," 1.00e15 1.00e17"}' load_segs.dat > load.dat
	rm one.dat two.dat 
# 
#   make a load for the trench model and force the grid into the gravity coordinate system
#	-R132/152/26.9923155926/33.0066161363
 	rm fz_xy.grd
 	vload 263.0 -40.0039553203 load.dat fz_xy.grd
    cp fz_xy.grd fz_merc.grd
# 
#   keep the fz.grd file in the Mercator space and the 
#   fz_xy.grd in the xy space
# 
	grdedit fz_merc.grd `grdinfo -I- grav_merc.grd`
#  
#	grdmath fz_xy.grd 0.0 MUL -4.0e-6 ADD = curv0_xy.grd
	grdmath fz_xy.grd 0.0 MUL -2.0e-8 ADD = curv0_xy.grd 

	cp age_merc.grd age_xy.grd
	grdedit age_xy.grd `grdinfo -I- fz_xy.grd`
# 
# 
#    APPLY MASK (assume xy file for mask polygon exists, can convert from KML later)
#