[cig-commits] r14015 - in mc/2D/ConMan/trunk/visual: . gmt

[becker at geodynamics.org]

Author: becker
Date: 2009-02-03 17:23:29 -0800 (Tue, 03 Feb 2009)
New Revision: 14015

Added:
   mc/2D/ConMan/trunk/visual/gmt/
   mc/2D/ConMan/trunk/visual/gmt/split_field_out
Log:
Added simple GMT plotting script.




Added: mc/2D/ConMan/trunk/visual/gmt/split_field_out
===================================================================
--- mc/2D/ConMan/trunk/visual/gmt/split_field_out	                        (rev 0)
+++ mc/2D/ConMan/trunk/visual/gmt/split_field_out	2009-02-04 01:23:29 UTC (rev 14015)
@@ -0,0 +1,89 @@
+#!/bin/bash
+#
+# splits conman ascii output into GMT/netcdf grd files
+#
+model=${1-new}			# model name
+plot=${2-3}			# 0: produce grids 1: grids and plot 2: plot and remove grids 3: convert to GIF
+
+
+if=field.$model			# input file
+op=grd.$model			# output prefix
+
+
+tmpn=/tmp/$USER.$HOST.$$.torque_tmp # temp storage
+trap "rm -f $tmpn.* ; exit" 0 1 2 15
+
+
+pname=`basename $0`
+
+
+if [ ! -s $if ];then
+    echo $pname: cannot find conman field output $if
+    exit
+fi
+
+head -1 $if > $tmpn.dat
+read nsd nx nz np nstep time < $tmpn.dat
+
+((nline=np+2))
+
+split -d --suffix-length=5 --lines $nline $if $tmpn.t.
+
+total_step=`ls $tmpn.t.* | wc -l`
+echo $pname: determined $nstep timesteps in model $model
+oc=0
+while [ $oc -lt $total_step ];do
+    tlabel=`echo $oc | gawk '{printf("%05i",$1)}'` # timelabel
+    fname=$tmpn.t.$tlabel
+    #echo $pname: working on $fname
+    if [ ! -s $fname ];then
+	echo $pname: output error, did not find split file $fname
+	exit
+    fi
+    head -1 $fname > $tmpn.dat
+    read nsd nelx nely np nstep time < $tmpn.dat # read time
+    if [ $oc -eq 0 ];then	# determine geometry for first split file
+	minmax -H2 -C $fname > $tmpn.dat
+	read nmin nmax xmin xmax ymin ymax vxmin vxmax vymin vymax tmin tmax < $tmpn.dat
+	reg=-R$xmin/$xmax/$ymin/$ymax
+	inc=`echo $xmin $xmax $nelx $ymin $ymax $nely | gawk '{printf("-I%.6f/%.6f",($2-$1)/($3),($5-$4)/($6))}'`
+    fi
+    # extract grids
+    tc=4
+    for t in vx vy t;do
+	ofile=$op.$t.$tlabel
+	gawk -v tc=$tc '{if(NR>2)print($2,$3,$(tc))}' $fname | \
+	    xyz2grd -G$ofile  $reg $inc -Dremark="t = $time"
+	((tc=tc+1))
+    done
+    if [ $plot -gt 1 ];then
+	if [ $oc -eq 0 ];then
+	    makecpt -Cpolar -T0/1/.1 > $tmpn.cpt
+	    proj=-Jx3
+	fi
+
+	ofile=$op.$tlabel.ps
+	grdimage -C$tmpn.cpt $op.t.$tlabel -P $proj $reg -Ba.5f.05:"x":/a.5f.05:"z":WesN -K > $ofile
+	grdvector $op.vx.$tlabel  $op.vy.$tlabel  $proj $reg -Q0.015/0.06/0.045 -N \
+	    -I0.2 -S2000 -O -K -Ggray -W0.5 >> $ofile
+	psscale -N50 -D0.2/-.2/2/.2h -C$tmpn.cpt -B.25/:"T": -O  >> $ofile
+	echo $pname: written to $ofile
+	if [ $plot -eq 3 ];then	# convert
+	    echo $pname: converting to $op.$tlabel.gif, deleting PS
+	    /usr/bin/convert $ofile $op.$tlabel.gif
+	    rm $ofile
+	fi
+    fi
+    if [ $plot -ge 2 ];then
+	rm $op.t.$tlabel $op.vx.$tlabel $op.vy.$tlabel 
+    fi
+    ((oc=oc+1))
+done
+
+if [ $plot -eq 3 ];then		# combine to movie
+    gifsicle $op.*.gif $model.movie.gif
+    rm $op.*.gif
+    echo $pname: output in $model.movie.gif
+
+
+fi