:

# Runs seasonal analysis on an Aeolus output data set
# $1 aeolus model output
# $2 FMS atmos model down output
# $3 FMS flux output

mydir="`dirname $0`"

usage() {
    echo "usage:"
    echo "$0 [aeolus model output] [atmos model down] [flux]"
    echo "Runs seasonal analysis on an Aeolus output data set"
    exit 1
}

usepyf=""
nodisplay=""
program="`which pyferret 2>/dev/null`"
if [ "X" = "X$program" ]
then program="`which ferret`"
else
    usepyf=1
    nodisplay=-nodisplay
fi


yseasmean() {
    # Prepare seasonal mean of a data file
    # $1 input data file name
    # $2 output file name

    # Set climatology taxis, such that Aeolus and Era have same
    # time axis for the climatologies and can be compared.
    # Also set vertical axis to pressure axis to enable comparison
    # with CM2 output and ERA-Interim.

    # First figure out number of levels, to select the apropriate
    # z-axis description file.

    setzaxis=""
    nlevels=`ncdump -h $1 | awk '/level = / { print $3;}'`
    if [ -f $mydir/zaxis_pressure_$nlevels ]
    then
	aeoluslevels=$nlevels
	setzaxis="setzaxis,$mydir/zaxis_pressure_$nlevels"
    else
	echo $1 has $nlevels vertical levels but there is no apropriate $mydir/zaxis_pressure_$nlevels
    fi

    cdo $setzaxis \
        -settaxis\,2013-02-01\,00:00\,3months \
	-yseasmean \
	-seldate\,0052-01-01T00:00:00\,0104-12-31T23:59:59 \
	$1 $2
    #ncatted -a axis\,level\,c\,c\,Z $2
    # cdo has added level bounds as 2-D variable,
    # but ferret needs a 1-D dimension variable. Sigh.
    # Copying from the original file is useless, because those
    # are still in meters, not pressure units.
    #ncks -A -v levelb $1 $2
    #ncatted -a edges\,level\,c\,c\,levelb $2
}

if [ ! -f "$3" ] ; then usage ; fi

tmpaeolus=Run_Seasonal-Aeolus-tmp.nc
tmpdown=Run_Seasonal-atmosdown-tmp.nc
tmpflux=Run_Seasonal-flux-tmp.nc

yseasmean "$1" $tmpaeolus
yseasmean "$2" $tmpdown
yseasmean "$3" $tmpflux

#exit

FER_GO="$FER_GO $mydir"

script=script-$$.jnl
tstamp="`date '+%b %Y'`"

echo 'can mode logo' > $script
#aeoluslevels=`ncdump -h $tmpaeolus | awk '/plev = / { print $3;}'`
aeoluslevels=6 # hardcoded for now
case "$aeoluslevels" in
6)
echo 'DEF AXIS/Z/UNITS=mbar/BOUNDS/DEPTH zax = {24,325,417,535,687,882}, {0,287, 287,368, 368,472, 472,607, 607,779, 779,1000}' >> $script
;;
17)
echo 'DEF AXIS/Z/UNITS=mbar/BOUNDS/DEPTH zax = {16,144,163,185,210,238,269,305,346,392,444,503,570,646,732,829,939}, {0,135, 135,153, 153,174, 174,197, 197,223, 223,253, 253,287, 287,325, 325,368, 368,417, 417,472, 472,535, 535,607, 607,687, 687,779, 779,882, 882,1000}' >> $script
;;
*) echo "Dont know how to create Z AXIS DEFINITION for $aeoluslevels levels"
   exit 42
   ;;
esac
echo 'go Albedo.jnl 1,DJF,"'$tmpaeolus'","'$tmpdown'","'$tmpflux'"' >> $script
echo 'go Albedo.jnl 2,MAM,"'$tmpaeolus'","'$tmpdown'","'$tmpflux'"' >> $script
echo 'go Albedo.jnl 3,JJA,"'$tmpaeolus'","'$tmpdown'","'$tmpflux'"' >> $script
echo 'go Albedo.jnl 4,SON,"'$tmpaeolus'","'$tmpdown'","'$tmpflux'"' >> $script
echo 'go Albedo.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'","'$tmpdown'","'$tmpflux'"' >> $script
echo 'go Energytransport.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Energytransport.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Energytransport.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Energytransport.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Energytransport.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
#echo quit >> $script
#ferret -gif -nojnl -script $script
#exit 0
echo 'go Heatflux.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Heatflux.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Heatflux.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Heatflux.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Heatflux.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go LWR_Radiation.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go LWR_Radiation.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go LWR_Radiation.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go LWR_Radiation.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go LWR_Radiation.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go Temperature.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Temperature.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Temperature.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Temperature.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Temperature.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go Watertransport.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Watertransport.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Watertransport.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Watertransport.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Watertransport.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go Clouds.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Clouds.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Clouds.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Clouds.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Clouds.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
#echo '!go Clouds2.jnl 1,DJF,"'$tmpaeolus'"' >> $script
#echo '!go Clouds2.jnl 2,MAM,"'$tmpaeolus'"' >> $script
#echo '!go Clouds2.jnl 3,JJA,"'$tmpaeolus'"' >> $script
#echo '!go Clouds2.jnl 4,SON,"'$tmpaeolus'"' >> $script
#echo '!go Clouds2.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go Humidity.jnl 1,DJF,"'$tmpaeolus'","'$tmpflux'"' >> $script
echo 'go Humidity.jnl 2,MAM,"'$tmpaeolus'","'$tmpflux'"' >> $script
echo 'go Humidity.jnl 3,JJA,"'$tmpaeolus'","'$tmpflux'"' >> $script
echo 'go Humidity.jnl 4,SON,"'$tmpaeolus'","'$tmpflux'"' >> $script
echo 'go Humidity.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'","'$tmpflux'"' >> $script
#echo 'go SWR_Radiation.jnl 1,DJF,"'$tmpaeolus'","'$tmpdown'"' >> $script
#echo 'go SWR_Radiation.jnl 2,MAM,"'$tmpaeolus'","'$tmpdown'"' >> $script
#echo 'go SWR_Radiation.jnl 3,JJA,"'$tmpaeolus'","'$tmpdown'"' >> $script
#echo 'go SWR_Radiation.jnl 4,SON,"'$tmpaeolus'","'$tmpdown'"' >> $script
#echo 'go SWR_Radiation.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go SWR_Radiation2.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go SWR_Radiation2.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go SWR_Radiation2.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go SWR_Radiation2.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go SWR_Radiation2.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo 'go Wind.jnl 1,DJF,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind.jnl 2,MAM,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind.jnl 3,JJA,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind.jnl 4,SON,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind_500mb.jnl 1,DJF,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind_500mb.jnl 2,MAM,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind_500mb.jnl 3,JJA,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind_500mb.jnl 4,SON,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Wind_500mb.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'","'$tmpdown'"' >> $script
echo 'go Precipitation.jnl 1,DJF,"'$tmpaeolus'"' >> $script
echo 'go Precipitation.jnl 2,MAM,"'$tmpaeolus'"' >> $script
echo 'go Precipitation.jnl 3,JJA,"'$tmpaeolus'"' >> $script
echo 'go Precipitation.jnl 4,SON,"'$tmpaeolus'"' >> $script
echo 'go Precipitation.jnl 1:4@ave,ANNUAL,"'$tmpaeolus'"' >> $script
echo quit >> $script
#ferret -batch -nojnl -script $script
#*** NOTE: Only one window can be open in batch mode
$program $nodisplay -nojnl -script $script

#rm -f $script $tmpaeolus $tmpdown $tmpflux *.plt*