#!/bin/sh
# 2006_0917 NISHIDA, Kenlo
# run this under GRASS
# i=2006_05_molpmonth
gmtset LABEL_FONT_SIZE      18
gmtset ANNOT_FONT_SIZE_PRIMARY 18
gmtset ANNOT_FONT_SIZE_SECONDARY 18
gmtset DOTS_PR_INCH        300

# map region. "r" is the resolution (in degree)
r=`g.region -g | grep "^nsres=" | awk 'BEGIN{FS="="}{print $NF}'`
e=`g.region -g | grep "^e=" | awk 'BEGIN{FS="="}{print $NF}' r=$r`
w=`g.region -g | grep "^w=" | awk 'BEGIN{FS="="}{print $NF}' r=$r`
s=`g.region -g | grep "^s=" | awk 'BEGIN{FS="="}{print $NF}' r=$r`
n=`g.region -g | grep "^n=" | awk 'BEGIN{FS="="}{print $NF}' r=$r`
scale=0.40
yymm=`echo $i | sed 's/_/ /g' | awk '{print $1"/"$2}'`

# export map data from GRASS. Anomaly is converted into "percent" scale.
r.out.bin input=$i output=$i.raw
r.mapcalc "dummy=${i}_anom*100.0"
r.out.bin input=dummy output=${i}_anom.raw
g.remove rast=dummy

# convert to GMT-grd format
xyz2grd $i.raw -Dm/m/m/1/0/=/= -G${i}.grd -R${w}/${e}/${s}/${n} -I${r} -V -ZTLd -F
xyz2grd ${i}_anom.raw -Dm/m/m/1/0/=/= -G${i}_anom.grd -R${w}/${e}/${s}/${n} -I${r} -V -ZTLd -F

### figure of mol/month
# color table for mol/month
cat <<EOF > molpmonth.cpt
-10000   0   0   0            0      0   0   0
0        0   0  70          400      0   0 255
400      0   0 255          800      0 255   0
800      0 255   0         1200    255 255   0
1200   255 255   0         1600    255  70  70
1600   255  70  70         2000     70   0   0
2000     0   0   0        10000      0   0   0
EOF

out=${i}.ps
psbasemap -Jm${scale} -R${w}/${e}/${s}/${n} -B10f5 -X4 -Y7 -P -K -V -U$0  > $out
grdimage ${i}.grd -R -Jm${scale} -Cmolpmonth.cpt -O -P -K -V>> $out
#pscoast -Jm${scale} -R -Lf136.5/15/33/1000 -Dh -N1/2/5/5/5 -W2/5/5/5 -P -O -K -V >> $out
pscoast -Jm${SCALE} -R -Dh -P -V -K -N1/1/50/50/50 -W2/20/20/20 -O >> $out
echo "131 28 20 0 0 0 Monthly PAR "$yymm | pstext -Jm -R -O -K >> $out
echo "145 39 20 0 0 0 mol" | pstext -Jm -R -O -K >> $out
cat molpmonth.cpt | awk '0<=$1 && $1<2000 {print $0}' > dummy.cpt
psscale -Cdummy.cpt -D8.5c/3.6c/3.5c/0.3c -O >> $out
ps2pdf $out
convert $out -rotate 90 ${i}.png

### figure of anomaly
# color table for anomaly (percent)
cat <<EOF > anomaly.cpt
-50     80  80 180          -25     80  80 255
-25     80  80 255            0    255 255 255
  0    255 255 255           25    255  80  80
 25    255  80  80           50    180  80  80
EOF
out=${i}_anom.ps

psbasemap -Jm${scale} -R${w}/${e}/${s}/${n} -B10f5 -X4 -Y7 -P -K -V -U$0  > $out
grdimage ${i}_anom.grd -R -Jm${scale} -Canomaly.cpt -O -P -K -V>> $out
#pscoast -Jm${scale} -R -Lf136.5/15/33/1000 -Dh -N1/2/5/5/5 -W2/5/5/5 -P -O -K -V >> $out
pscoast -Jm${SCALE} -R -Dh -P -V -K -N1/1/50/50/50 -W2/20/20/20 -O >> $out
echo "131 28 20 0 0 0 PAR anomaly "$yymm | pstext -Jm -R -O -K >> $out
echo "145 39 20 0 0 0 %" | pstext -Jm -R -O -K >> $out
psscale -Canomaly.cpt -D8.5c/3.6c/3.5c/0.3c -O >> $out
ps2pdf $out
convert $out -rotate 90 ${i}_anom.png

convert ${i}.png ${i}_anom.png -append -rotate 270  ${i}_2.png