#!/bin/bash
#
# given a spherical harmonics model in layer format, create a gridded representation
# uses tools from the shansyn package, 
# http://www-udc.ig.utexas.edu/external/becker/sdata.html
#
# $Id: sh_to_grds,v 1.3 2016/07/27 02:28:30 becker Exp $
#
mfile=$1	# SH model file, e.g. savani.31.m.ab, required
dx=${2-0.5}	# grid spacing for the global grids, with default
gname=${3-dv}	# prefix for the grid file names, dv.1.grd dv.2.grd etc.
dfile=${4-depths.dat}		# file with depths, as given in m.ab file, bottom up
#
#
#
if [ $# -lt 1 ];then
    echo $0: usage
    echo
    echo `basename $0` model.lmax.m.ab \[dx, $dx\] \[gname, $gname\] \[dfile, $dfile\] 
    echo
    echo converts SH model file model.lmax.m.ab into GMT grds of dx grid spacing
    echo named grdname.\*.grd and prints depth to dfile
    echo
    exit
fi
if [ ! -s $mfile ];then
    echo $0: model file $mfile not found
    exit
fi

ncpus=${NR_CPUS-1}		# how many CPUs to use 

tmpn=/tmp/$USER.$HOST.$$.stg	# temp files
trap "rm -f $tmpn.* ; exit" 0 1 2  15
extract_model_depths $mfile | gawk '{print($2)}' > $dfile # depth file
depths=`gawk '{printf("%s ",$1)}' $dfile`
n=`wc -l $dfile`
echo $0: written $n layer depths to $dfile, expanding at $dx spacing
j=0;i=1
for d in $depths;do
    echo $0: expanding layer $i at depth $d
    (extract_layer $mfile $d | shsyn $dx 0 $gname.$i ) &
    ((j=j+1))
    if [ $j -eq $ncpus ];then	
	wait
	j=0
    fi
    ((i=i+1))
done
wait
echo $0: written to $gname.*.grd
rm expansion.par
