#!/bin/bash
#
D=../01_GPS/wus_gps_final_crowell.dat
R=-R124.5W/114.0W/31.0N/40.5N
#
#  clean out old files
#
rm -f data.lluv junk* tmp*
#
# reformat the data and set the minimum sigma at 0.6 mm/yr
#
gmt select $D $R -fg | awk '{su = ($5 < .6 ? .6 : $5);sv = ($6 < .6 ? .6 : $6);print($1,$2,$3,$4,su,sv) }' > data.lluv
#
#  extract the model from the grid
#
gmt grdtrack data.lluv -G../02_models/$1/ve_tot.grd > junk1.dat
gmt grdtrack junk1.dat -G../02_models/$1/vn_tot.grd > junk2.dat
#
# gridding the clean.dat
#
awk '{print ($1,$2,$3-$7)}' < junk2.dat | gmt blockmedian $R -I1m -r | awk '{print ($1, $2, $3)}'  >  GPS_resid_east.dat
awk '{print ($1,$2,$4-$8)}' < junk2.dat | gmt blockmedian $R -I1m -r | awk '{print ($1, $2, $3)}'  >  GPS_resid_north.dat
#
# combine the east and north residuals
#
gmt convert -A GPS_resid_east.dat GPS_resid_north.dat -o0-2,5 > blk.lluv
#
#  do the gridding
#
gmt gpsgridder `gmt grdinfo -I ../02_models/$1/ve_tot.grd` `gmt grdinfo -I- ../02_models/$1/ve_tot.grd` -Gtmp_%s.nc blk.lluv -fg -Vl -Fd6 -r -Cn1800 -S0.50
#
#  now grid the residuals add back to the velocity and resample to the same grid
#
gmt grdmath ../02_models/$1/ve_tot.grd tmp_u.nc ADD = tmp_e.grd
gmt grdmath ../02_models/$1/vn_tot.grd tmp_v.nc ADD = tmp_n.grd 
grdmath topo.grd -100 GT 0. NAN = tmp1.grd
grdsample tmp1.grd -R-123/-114/32.0/37.5 -I0.01 -Gmask.grd
grdsample tmp_e.grd -R-123/-114/32.0/37.5 -I0.01 -Gtmp_e2.grd
grdsample tmp_n.grd -R-123/-114/32.0/37.5 -I0.01 -Gtmp_n2.grd
grdmath tmp_e2.grd mask.grd OR = ve_$1.grd
grdmath tmp_n2.grd mask.grd OR = vn_$1.grd

./plot_Lat_Long.com -123 -114 32.0 37.5 2.00 vn_$1 vn 5
./plot_Lat_Long.com -123 -114 32.0 37.5 2.00 ve_$1 ve 5
#
rm GPS_resid_east.dat GPS_resid_north.dat
rm junk*
#
# first take the gradidient of the velocity grids
#
gmt grdgradient ve_$1.grd -Gvxx.grd -A270 -fg -V
gmt grdgradient ve_$1.grd -Gvxy.grd -A180 -fg -V
gmt grdgradient vn_$1.grd -Gvyx.grd -A270 -fg -V
gmt grdgradient vn_$1.grd -Gvyy.grd -A180 -fg -V
#
gmt grdmath vxx.grd 1000. DIV 1.e9 MUL = exx_$1.grd
gmt grdmath vyy.grd 1000. DIV 1.e9 MUL = eyy_$1.grd
gmt grdmath vxy.grd vyx.grd ADD 2000. DIV 1.e9 MUL = exy_$1.grd
#