# script to project model velocity files (vx, vy, vz) from PoD to geographic coordinates


# crop the grids

grdcut xvel2011.grd -Gxjunk.grd -R0/1000/0/1700 -V
grdcut yvel2011.grd -Gyjunk.grd -R0/1000/0/1700 -V
grdcut zvel2011.grd -Gzjunk.grd -R0/1000/0/1700 -V

# grd2xyz for an ascii table

grd2xyz xjunk.grd > xv.xyz
grd2xyz yjunk.grd > yv.xyz
grd2xyz zjunk.grd > zv.xyz

# backshift x- ( -400 )  and y- (-850) km  from model 0/0 center
awk '{print ($1-400), ($2-850), $3}' < xv.xyz > xvs.xyz
awk '{print ($1-400), ($2-850), $3}' < yv.xyz > yvs.xyz
awk '{print ($1-400), ($2-850), $3}' < zv.xyz > zvs.xyz


# prepare the file for trans_pole

# print ykm xkm vy vx 00 00 00 for trans_pole velocity (vx, vy)
# print ykm xkm vz for trans_pole scalar (vz)

paste xvs.xyz yvs.xyz zvs.xyz > xyzvs.xyz
awk '{print $2, $1, $6, $3, 0.00, 0.00, 0.00}' < xyzvs.xyz > vxvy.xyz
awk '{print $2, $1, $9}' < xyzvs.xyz > vz.xyz


# transform back to geographic coordinates

trans_pole -2 2 285.6 50.1 276 54 < vxvy.xyz > lonlatvevn.dat
trans_pole -2 1 285.6 50.1 276 54 < vz.xyz > lonlatvu.dat

# should give Podlong Podlat vy vx 0 0 0 lon lat ve vn se sn 0 
# should give Podlong Podlat vz lon lat  

awk '{print -(360-$8), $9, $10}' < lonlatvevn.dat > ve.xyz
awk '{print -(360-$8), $9, $11}' < lonlatvevn.dat > vn.xyz
awk '{print -(360-$4), $5, $3}' < lonlatvu.dat > vu.xyz


# next run gmt surface program for x and y surface interpolation

surface ve.xyz -Gvjunk.grd -R-127/-110.3330/28/41.5 -I.01/.01 -V
grdsample vjunk.grd -Gve.grd -T

surface vn.xyz -Gvjunk.grd -R-127/-110.3330/28/41.5  -I.01/.01 -V
grdsample vjunk.grd -Gvn.grd -T

surface vu.xyz -Gvjunk.grd -R-127/-110.3330/28/41.5 -I.01/.01 -V
grdsample vjunk.grd -Gvu.grd -T

 ./plot_Lat_Long.com -125 -114 30.0 41.5 0.63 vn 5
 ./plot_Lat_Long.com -125 -114 30.0 41.5 0.63 ve 5 
 ./plot_Lat_Long.com -125 -114 30.0 41.5 0.63 vu 1