#!/bin/bash
#
# given a gridded model in GMT/Netcdf grd layers, create a spherical harmonics expansion model
#
# uses tools from the shansyn package, 
# http://www-udc.ig.utexas.edu/external/becker/sdata.html
#
# $Id: grds_to_sh,v 1.2 2016/07/27 02:36:40 becker Exp $
#
gname=$1     	# prefix for the grid file names, dv.1.grd dv.2.grd etc
lmax=${2-31}	# max spherical harmonics degree
dfile=${3-depths.dat}		# file with depths, bottom up
mfile=${4-$gname.$lmax.m.ab}	# SH model output file name 
#
#
ncpus=${NR_CPUS-1}		# how many CPUs to use 
#
#
#
if [ $# -lt 1 ];then
    echo $0: usage
    echo
    echo `basename $0` gname \[lmax, $lmax\] \[dfile, $dfile\] \[mfile, gname.$lmax.m.ab\] 
    echo
    echo converts gname.1.grd gname.2.grd ... grids to a spherical harmonics
    echo model file mfile, using max degree lmax, and expecting depth given bottom up
    echo in dfile
    echo
    exit
fi
tmpn=/tmp/$USER.$HOST.$$.stg	# temp files
trap "rm -f $tmpn.* ; exit" 0 1 2  15


if [ ! -s $dfile ];then
    echo $0: depth file $dfile not found
    exit
fi
depths=`gawk '{printf("%s ",$1)}' $dfile`
n=`echo $depths | gawk '{print(NF)}'`
echo $0: expanding $gname.\?.grd to lmax $lmax
echo $0: read $n depths from $dfile

# expand
i=1;j=0
for d in $depths;do
    if [ ! -s $gname.$i.grd ];then
	echo $0: grid $gname.$i.grd not found
	exit
    fi
    echo $0: working on $gname.$i.grd at depth $d
    ( shana $lmax $gname.$i  > $tmpn.$i.ab ) &
    ((j=j+1))
    if [ $j -eq $ncpus ];then	
	wait
	j=0
    fi
    ((i=i+1))
done
wait
# assemble
echo $n > $mfile
i=1;j=0
for d in $depths;do
    echo $d >> $mfile
    if [ ! -s $tmpn.$i.ab ];then
	echo $0: error, expansion for depth $d index $i not produced
	exit
    fi
    cat $tmpn.$i.ab >> $mfile
    ((i=i+1))
done
echo $0: written to $mfile


