#!/bin/csh
# cut and plot the model
# 

set par = $1 		## Parameter to be plotted (S/X)
set dir = /data/08/federica/WINV/     ## model directory
set nmodel  = $2	## Model name
set A2index = $3	## Level

set ifrmean = 1		## 1/0 Do/don't remove mean value
set step = 2.5		## Grid size in degree for extracting and plotting

set basef = /data/08/federica/ACAL_SPL/grid/grid.$A2index

##############################################

PLOT:

set psfile = $nmodel\_$par\_na.ps

\rm  lu_*

if ( ! ( -f $dir/$nmodel ) ) then
ll $dir/$nmodel
exit
endif

awk '{if(NR!=1)print$0}' $basef > tmp.grid
ln -s $basef lu_hknots
ln -s $dir/$nmodel lu_coef

if ( $par != 'M') then
## change depth  HERE
#Source code: /data/08/federica/SRC/SPH/ANISO/grid_mA3d_anis_na.c
/data/08/federica/bin/grid_mA3d_anis_na.xx $step <<END
$par
75
100
150
200
250
300
350
400
-1
END
endif

if ($par == 'M') then
/data/08/federica/bin/grid_mA3d_anis_na.xx $step <<END
$par
END
endif
echo grid_A3d_anis.xx DONE

if ( $par != 'M') then
set xyzout = grid_A3d_anis.out 
set depinfo = `awk '{if(NR==1)print $0}' $xyzout `
foreach d ( $depinfo )
echo $d | awk '{printf("%d \t",$1)}' >> tmp.$$
end
set deps = `cat tmp.$$ `
\rm -rf tmp.$$
 
echo DEPTH = $deps

set i = 0
foreach dep ( $deps )
@ i ++
@ nr = 2 + $i
awk '{if(NR>1) print $1,$2,$ND}' ND=$nr $xyzout > A3.$dep.xyz
echo A3.$dep.xyz ready 
end

if ( $ifrmean ) then
foreach dep ( $deps )
/data/07/gung/,sphA1/,test/rmean.exe <<END
A3.$dep.xyz
A3.$dep.rmean.xyz
END
end
endif

#################### GMT #######################

gmtset BASEMAP_TYPE PLAIN
gmtset MEASURE_UNIT   cm
gmtset PAPER_MEDIA Letter
gmtset FRAME_PEN	 0.5p
gmtset FRAME_WIDTH	 0.1c
gmtset LABEL_FONT_SIZE	 10p
gmtset ANOT_FONT_SIZE	 10p
gmtset TICK_LENGTH	 0.05i

set xl = 9
set yl = 4.5
set xo = 3.5
set yo = 21.0
set mid = 180

set r = -R-130/20/-45/60r
#set r = -R-112/23/-70/45r
set projq = -JQ$mid/$xl
set proj = -JS-100/20/5.

set projx       = -JX$xl/$yl
set rx          = -R0/360/-90/90

set lsca        = -2/-1/6/0.4h

set grid        = swne
set area	= 6000/1
set cptf        = green4.cpt

set plate = /data/07/gung/,MODELS/plate.pacific.asc
set geology = /a/perry/data/08/federica/INPUT/GMT/na_geology
#set amps = ( 4 4 3.5 3 2.5 2 2 2 )
set amps = 4
/data/08/federica/SCRIPTS/make_green7z.s $amps > $cptf 
#set cptf = /data/08/federica/INPUT/COLOR_SCALE/color.craton.profile

@ i = 0
foreach dep ( $deps )
@ i ++

set f = A3.$dep.xyz
if ( $ifrmean ) then
set f = A3.$dep.rmean.xyz
endif

set info = `minmax -C $f | awk '{printf("%.2f %.2f %.2f\n", (-$5+$6)/2.*0.8, $5,$6)}'`

set step1 = 1.0

#if ( $dep == 150 && $par == 'X' ) then
#awk '{if ($1>180) $1=$1-360; print $1,$2,(1.0337*(1+($3/100.))-1)*100.}' $f > temp.dat
#else if ( $dep == 100 && $par == 'X' ) then
#awk '{if ($1>180) $1=$1-360; print $1,$2,(1.06*(1+($3/100.))-1)*100.}' $f > temp.dat
#else if ( $dep == 75 && $par == 'X' ) then
#awk '{if ($1>180) $1=$1-360; print $1,$2,(1.0733*(1+($3/100.))-1)*100.}' $f > temp.dat
#else if ( $dep == 200 && $par == 'X' ) then
#awk '{if ($1>180) $1=$1-360; print $1,$2,(1.008*(1+($3/100.))-1)*100.}' $f > temp.dat
#else
awk '{if ($1>180) $1=$1-360; print $0}' $f > temp.dat
#endif

mv temp.dat $f
surface $f -Gtmp.grd -I$step1/$step1 -T0.3 -R-180/0/10/90
#xyz2grd $f -Gtmp.grd -I$step1/$step1 -R-180/0/10/90 -V

if ( $i == 1 )  then
grdimage tmp.grd $proj $r -X$xo -Y$yo -C$cptf -K -P >! $psfile
psxy -: $geology -R -JS -W7 -M -K -O -V >> $psfile
psxy tmp.grid -R -JS -Sc0.1 -G0/0/0 -K -O -V >> $psfile
pstext <<END $projq $rx -K -N -O >>$psfile
0 130 10 0 0 1 $nmodel 
END
else if ( $i == 5 ) then
grdimage tmp.grd $proj $r -Y16.2 -X9.0 -C$cptf -K -O -P >> $psfile
else
grdimage tmp.grd $proj $r -Y-5.4 -C$cptf -K -O -P >> $psfile
endif 

pstext $projq $rx -K -N -O  <<END>> $psfile
100 95 10 0 0 2 $dep km ( $info[2]/$info[3] )
END

pscoast -B$grid $r $proj -Dc -A$area -Na -W3/0 -O -K >> $psfile

psxy -: $plate $r $proj -M -W2t10:10 -O -K >> $psfile

end

if ( $par == 'S') then
#psscale -C$cptf -D$lsca -Ba2:"dlnv@-S@-(%)": -O -L >> $psfile
psscale -C$cptf -D$lsca -Ba2:"dlnv@-S@-(%)": -O >> $psfile
else
#psscale -C$cptf -D$lsca -Ba2:"dln@~x@~(%)": -O -L >> $psfile
psscale -C$cptf -D$lsca -Ba2:"dln@~x@~(%)": -O >> $psfile
endif

#\rm A3.*.xyz
\rm temp tmp.grd tmp.grid lu* bu* green4.cpt grid_A3d_anis.out
\rm temp tmp.grd tmp.grid lu* bu* grid_A3d_anis.out
endif

if ( $par == 'M') then
set solution = grid_A3d_anis.out 
awk '{if ($1>180) $1=$1-360; print $0}' $solution > temp.dat
mv temp.dat $solution

gmtset COLOR_NAN 190/190/190
gmtset FRAME_PEN 1.25
gmtset ANOT_FONT_SIZE 13p
gmtset LABEL_FONT_SIZE 13p
gmtset BASEMAP_TYPE PLAIN
gmtset MEASURE_UNIT inch
gmtset PAPER_MEDIA Letter

set xl = 9
set yl = 4.5
set xo = 1.2
set yo = 3.0
set mid = 180

set r = -R-130/20/-45/60r
set projq = -JQ$mid/$xl
set proj = -JS-100/20/6.

set projx       = -JX$xl/$yl
set rx          = -R0/360/-90/90

set lsca        = 3/-0.7/5/0.3h

set grid        = swne
set area	= 6000/1
set cptf        = /data/08/federica/INPUT/COLOR_SCALE/color.na

set plate = /data/07/gung/,MODELS/plate.pacific.asc
set step1 = 1.0
set info = `minmax -C $solution | awk '{printf("%.2f %.2f %.2f\n", (-$5+$6)/2.*0.8, $5,$6)}'`

surface $solution -Gtmp.grd -I$step1/$step1 -T0.3 -R-180/0/10/90
grdimage tmp.grd $proj $r -X$xo -Y$yo -C$cptf -K -P >! $psfile
pscoast -B$grid $r $proj -Dc -A$area -W3/0 -O -K >> $psfile
psxy -: $plate $r $proj -M -W2t10:10 -O -K >> $psfile

pstext $projq $rx -K -N -O  <<END>> $psfile
115 120 12 0 0 2 ( $info[2] / $info[3] )
END

psscale -C$cptf -D$lsca -Ba5:"Moho depth (km)": -O -L >> $psfile
\rm tmp.grd tmp.grid lu* bu* grid_A3d_anis.out
endif

echo $psfile

ghostview $psfile &