#!/bin/csh -x
# AWW-2006
# script to plot routed flows for observed & up to 3 simulations
# generalized max/intervals on y-axis
# 3 plots:  
#	1) monthly flow timeseries
#	2) monthly average for period

#input flow files have format: year month obs_data
#gmtset TICK_LENGTH 0.03

set LOC = "$1"
set OBSFLOW = "$2"
set RUN = "$3"

#set GMTPATH = "/opt/gmt-4.3.1/bin" # This is local to Mercury
#setenv LD_LIBRARY_PATH "/home1x/cics/pcic/opt/netcdf-3.6.3/lib" # And this.

set OBSFL = $OBSFLOW # /home1x/cics/pcic/vic/autocal/models/.BASIN_TEMPLATE/obs/$LOC.obs.month.kaf  # obs in kaf
set SIMFL = ../output/combined_flow_sacramento_tst/$LOC.month  # in cfs
#set SIMFL = $LOC.month
set BASEFL = /home/disk/columbia2/marketa/NPR/2d_simulations/output/base_4.1.2/routed_flows/$LOC.month

set OUTF = ./$LOC.month.ps
echo $OUTF

set TITLE = "$LOC ($RUN) net monthly mean inflow (cfs): OBS, OPTI, BASE"
set LBL = (OBS SIM_OPTI BASE)
set AVG_SYR = 1976  # can set averaging period if you want
set AVG_EYR = 1990  # need to change down below too

set YEARSTRT = 1976
set YEAREND = 1990

# ------ SET VARIABLES ------------------------------------------ 
set PROJ = -JX14/2.3
set COORD1 = -R1976/1979/0/
set COORD2 = -R1979/1981/0/
set COORD3 = -R1981/1983/0/
set COORD4 = -R1983/1985/0/
set COORD5 = -R1985/1988/0/
set COORD6 = -R1988/1990/0/
set tmp = `awk '{if($1>='$YEARSTRT' && $1<='$YEAREND') print $3}' $OBSFL | sort -n | tail -1`
set MAX1 = `/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/ascale 0 $tmp f l | awk '{print $2}' `
set INTV = `/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/ascale 0 $tmp f l | awk '{print $3}' `
echo timeseries tmp = $tmp MAX = $MAX1

echo -n pstext output lines:
GMT pstext -P -K -JX6.5/10 -R0/1/0/1 -N << END | tee $OUTF | wc -l
1.25 2.65 16 0 21 6 $TITLE
END

# -- print period 1 --------------------------------------------------

echo -n psxy output lines:
awk '{printf("%.3f %.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |   \
GMT psxy -O -K -B2f1/$INTV\:.:SWne $COORD1$MAX1 $PROJ -W3 -X1.3 -Y23 -V | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%.3f %.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |       \
GMT psxy -O -K $COORD1$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%.3f %.3f\n", $1+$2/12-1/24, $3)}' $BASEFL |       \
GMT psxy -O -K $COORD1$MAX1 $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

# -- print period 2 --------------------------------------------------

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |    \
GMT psxy -O -K -B2f1/$INTV\:.:SWne $COORD2$MAX1 $PROJ -W3 -Y-3.5 | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |\
GMT psxy -O -K $COORD2$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $BASEFL |\
GMT psxy -O -K $COORD2$MAX1 $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

# -- print period 3 --------------------------------------------------

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |    \
GMT psxy -O -K -B2f1/$INTV\:.:SWne $COORD3$MAX1 $PROJ -W3 -Y-3.5 | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |\
GMT psxy -O -K $COORD3$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $BASEFL |\
GMT psxy -O -K $COORD3$MAX1 $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

# -- print period 4 --------------------------------------------------

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |    \
GMT psxy -O -K -B2f1/$INTV\:.:SWne $COORD4$MAX1 $PROJ -W3 -Y-3.5 | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |\
GMT psxy -O -K $COORD4$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $BASEFL |\
GMT psxy -O -K $COORD4$MAX1 $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

# -- print period 5 --------------------------------------------------

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |    \
GMT psxy -O -K -B2f1/$INTV\:.:SWne $COORD5$MAX1 $PROJ -W3 -Y-3.5 | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |\
GMT psxy -O -K $COORD5$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l
echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $BASEFL |\
GMT psxy -O -K $COORD5$MAX1 $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

# -- print period 6 --------------------------------------------------

#echo -n psxy output lines:
#awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $OBSFL |    \
#   psxy -O -K -B2f1/$INTV\:.:SWne $COORD6$MAX1 $PROJ -W3 -Y-3.5 | tee -a $OUTF | wc -l
#echo -n psxy output lines:
#awk '{printf("%8.3f %8.3f\n", $1+$2/12-1/24, $3)}' $SIMFL |\
#   psxy -O -K $COORD6$MAX1 $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l

# %%%%%%%%%%%%%% now look at a mean monthly comparison %%%%%%%%%%%

# find time limits of simulation to get averages for
# uses later start date and earlier end date
rm -rf junk
head -1 $OBSFL > junk
head -1 $SIMFL >> junk
head -1 $BASEFL >> junk
set SYR = `sort -nr junk | awk '{if(NR==1)print $1}' `
rm -rf junk
tail -1 $OBSFL > junk
tail -1 $SIMFL >> junk
tail -1 $BASEFL >> junk
set EYR = `sort -n junk | awk '{if(NR==1)print $1}' `
echo $SYR $EYR summaries

rm -rf junk

# get monthly averages -------------------------------
/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/monTS_2_monAVG.scr $OBSFL $AVG_SYR $AVG_EYR -99.0 > obs.mon.avg
/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/monTS_2_monAVG.scr $SIMFL $AVG_SYR $AVG_EYR -99.0 > sim.mon.avg
/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/monTS_2_monAVG.scr $BASEFL $AVG_SYR $AVG_EYR -99.0 > base.mon.avg
cat obs.mon.avg sim*.mon.avg > junk

set tmp = `awk '{print $2}' junk | sort -n | tail -1`
set MAX2 = `/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/ascale 0 $tmp f l | awk '{print $2}' `
set INTV2 = `/home/disk/larix/marketa/NPR/autocal/mocom_npr_0.5deg_roagg/sources/ascale 0 $tmp f l | awk '{print $3}' `
set COORD_AVG = -R0.5/12.5/0/$MAX2

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1, $2)}' obs.mon.avg |    \
GMT psxy -O -K -B1:"month":/$INTV2\:.:SWne $COORD_AVG $PROJ -W3 -Y-3.85 | tee -a $OUTF | wc -l

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1, $2)}' sim.mon.avg |    \
GMT psxy -O -K $COORD_AVG $PROJ -W1/255/0/0:td | tee -a $OUTF | wc -l

echo -n psxy output lines:
awk '{printf("%8.3f %8.3f\n", $1, $2)}' base.mon.avg |    \
GMT psxy -O -K $COORD_AVG $PROJ -W1/0/0/255:td | tee -a $OUTF | wc -l

set HGT2 = `echo $MAX2 | awk '{print $1*1.25}' `
echo -n pstext output lines:
GMT pstext -O -K $COORD_AVG $PROJ -V -N << END | tee -a $OUTF | wc -l
#the first two numbers are all relative, so it must be set relative to the -B# set above 
6.5 $HGT2 14 0 21 6 Long Term Mean Daily Flow (CFS), $AVG_SYR-$AVG_EYR
END

# plot legend -----------------------------------------------------

echo -n pstext output lines:
GMT pstext -R0/1/0/1 $PROJ -N -K -O -Y-2.7 << END | tee -a $OUTF | wc -l
0.05 0.20 12 0 21 5 $LBL[1]
0.30 0.20 12 0 21 5 $LBL[2]
0.80 0.20 12 0 21 5 $LBL[3]
END

echo -n psxy output lines:
GMT psxy -O -K -N -R0/1/0/1 $PROJ -W3 << END | tee -a $OUTF | wc -l
0.00 0.20
0.03 0.20
0.06 0.20
END
echo -n psxy output lines:
GMT psxy -O -K  -N -R0/1/0/1 $PROJ -W1/255/0/0:td << END | tee -a $OUTF | wc -l
0.25 0.20
0.28 0.20
0.31 0.20
END
echo -n psxy output lines:
GMT psxy -O -N -R0/1/0/1 $PROJ -W1/0/0/255 << END | tee -a $OUTF | wc -l
0.75 0.20
0.78 0.20
0.81 0.20
END

echo PostScript file now contains `wc -l $OUTF` lines.  #DEBUG

#unsetenv LD_LIBRARY_PATH

\rm sim*.mon.avg obs.mon.avg base.mon.avg junk
#gv $OUTF


