#!/bin/sh
#=======================================================================
# Plot time averaged sst and sss anomalies               --- sst_sss ---
# $Id: sst_sss_jobdef 657 2012-02-22 23:48:49Z acrnrls $
#=======================================================================
#
# Adapted from deck sst_sss_avg_t63l33_from_mon_sst_sss_cent_on_pac.dk
# orig. from Bill Merryfield
#
# Larry Solheim May, 2007
#=======================================================================
# keyword :: sst_sss
# short description :: Plot time averaged sst and sss anomalies
set -a
. betapath2

#  * ........................... Parmsub Parameters ............................

username="acrnxxx"; user="XXX";
runid=job000; crawork="${runid}_job"

# These variables are set when the job string is created
run_start_year=NotSet
run_start_month=NotSet
run_stop_year=NotSet
run_stop_month=NotSet
current_year=NotSet
current_month=NotSet
previous_year=NotSet
previous_month=NotSet
next_year=NotSet
next_month=NotSet
if [ $current_year = "NotSet" ]; then
  echo "sst_sss: current_year is not set"
  exit 1
fi

# define $yearb and $yeare
if [ $previous_year = "NotSet" ]; then
  if [ $next_year = "NotSet" ]; then
    echo "sst_sss: Neither previous_year nor next_year are set"
    exit 1
  else
    yearb=$current_year   # memory99=1
    yeare=$next_year      # memory99=1
  fi
else
  yearb=$previous_year    # memory99=1
  yeare=$current_year     # memory99=1
fi

plunit="VIC"; plpfn=${runid}_${yearb}_$yeare
lopgm="lopgm_o4xa"
this_host=`hostname|cut -d'.' -f1`
case $this_host in
  # This is pollux, use 64 bit binaries
  ib3-*) lopgm=lopgm_o4xal ;;
esac

nextjob=on
noprint=on

jobname=sst_sss_; time="900" ; memory="900mb" ;

. comjcl.cdk

cat > Execute_Script <<'end_of_script'

  # bail is a simple error exit routine
  bail(){
    echo "sst_sss: "$1
    echo " " >> haltit
    exit 1
  }

  # files ${filem}_sst and ${filem}_sss must exist and contain monthly
  # averaged sst and sss from year $yearb month 1 to year $yeare month 12
  sst_sss_uxxx=sc;
  filem="${sst_sss_uxxx}_$runid"

  # lev is the level of the annually averaged sst and sss fields that
  # are found in the obs files.
  lev="   75"

  ccm="8   0";b="-2+";

  sstplt() {
    [ -z "$1" ] && bail "sstplt: missing start year"
    [ -z "$2" ] && bail "sstplt: missing end year"
    ptitle="ANN MEAN SST DIFFERENCE  ${runid} $1-$2 minus phc"
    ccc ggplot diff  <<EOF
   GGPLOT         -1   -1   -1$ccm      1.e0      -10.       10.     0.5000${b}8
$ptitle
          0  14  198  404  122  119  118  181  140
                 140  351  356  105  104  100  158                 
                -6.0      -5.0     -4.00      -3.0     -2.00      -1.0      0.00
                1.00       2.0      3.00       4.0      5.00       6.0     20.00
EOF
  }

  sssplt() {
    [ -z "$1" ] && bail "sssplt: missing start year"
    [ -z "$2" ] && bail "sssplt: missing end year"
    ptitle="ANN MEAN SSS DIFFERENCE  ${runid} $1-$2 minus phc"
    ccc ggplot diff  <<EOF
   GGPLOT         -1   -1   -1$ccm      1.e3      -10.        5.     0.2000${b}8
$ptitle
          0  14  198  404  122  119  118  181  140
                 140  351  356  105  104  100  158                 
                -3.0      -2.5     -2.00      -1.5     -1.00      -0.5      0.00
                0.50       1.0      1.50       2.0      2.50       3.0     10.00
EOF
  }

  # put TEMP and SALT from obs data into files TEMPobs and SALTobs

  # Ocean data description file
  datadesc="uwl_ogcm3.5_256x192l33v1_datadesc";
  access datadesc $datadesc
  grido datadesc betao dao dx dy dz
  release datadesc

  ccc select betao betaos << eor
SELECT    STEPS         0  88888888         1${lev}${lev}      BETA
eor

  # PHC 2.1 ann avg temp
  access datacr uwl_phc_p2.1_temp_ann_256x192l33v1
  rm -f datacs
  ccc select datacr datacs << eor
  SELECT. STEPS         0  99999999         1${lev}${lev}      TEMP
eor
  mlt datacs betaos TEMPobs
  release datacr

  # PHC 2.1 ann avg salinity
  access datacr uwl_phc_p2.1_salt_ann_256x192l33v1
  rm -f datacs
  ccc select datacr datacs << eor
  SELECT. STEPS         0  99999999         1${lev}${lev}      SALT
eor
  mlt datacs betaos SALTobs
  release datacr

  # access model time series data
  access TEMPmod ${filem}_sst
  access SALTmod ${filem}_sss

  # Put sst and sss anomalies averaged from yearb to yeare as
  # the first 2 plot frames in the file
  for VAR in sst sss; do
    case $VAR in
      sst) var="TEMP" ;;
      sss) var="SALT" ;;
      *) bail "Invalid variable $VAR" ;;
    esac

    # Add month and reformat for select
    y01=`echo $yearb|awk '{printf "%8d01",$1}' -`
    y12=`echo $yeare|awk '{printf "%8d12",$1}' -`

    # Use the level found in the first record with ibuf3=$var of
    # the model file as the level for the select done here.
    mlev=`ggstat ${var}mod|awk -v fld="$var" \
          '{if($4==fld){printf "%5d",$5;exit}}' -`

    # Extract the current year of monthly averages
    ccc select ${var}mod datamavgr <<EOF
  SELECT. STEPS${y01}${y12}         1${mlev}${mlev}      ${var}
EOF

    # Create average over the range yearb to yeare
    timavg datamavgr datamavgr1

    ggstat datamavgr1

    ccc relabl datamavgr1 datamavg  <<EOF
    RELABL GRID               
           GRID                                  0
EOF

    # determine difference from obs data
    sub datamavg ${var}obs diff

    ####### BUG Is this supposed to happen for sss
    if [ 1 -eq 2 -a "$VAR" = "sss" ] ; then
      echo "  XLIN        1000.0       0.0" > IN
      ccc xlin diff diffs input=IN
      echo "  XLIN           1.0       0.0" > IN
      ccc xlin diffs diff input=IN
    fi

    echo "$VAR difference for year $yyy"
    ggstat diff

    # Generate plot frame for ncgm file
    [ "$VAR" = "sst" ] && sstplt $yearb $yeare
    [ "$VAR" = "sss" ] && sssplt $yearb $yeare
  done

  if [ $yeare -ne $yearb ]; then
    # Add sst and sss annual average anomalies for each year
    # from yearb to yeare when yearb -ne yeare
    for VAR in sst sss; do
      case $VAR in
        sst) var="TEMP" ;;
        sss) var="SALT" ;;
        *) bail "Invalid variable $VAR" ;;
      esac
      # loop over years yearb to yeare
      yyy=`echo $yearb|awk '{printf "%3.3d",$1-1}' -`
      while [ $yyy -lt $yeare ]; do
        yyy=`echo $yyy|awk '{printf "%3.3d",$1+1}' -`

        # Add month and reformat for select
        y01=`echo $yyy|awk '{printf "%8d01",$1}' -`
        y12=`echo $yyy|awk '{printf "%8d12",$1}' -`

        # Use the level found in the first record with ibuf3=$var of
        # the model file as the level for the select done here.
        mlev=`ggstat ${var}mod|awk -v fld="$var" \
              '{if($4==fld){printf "%5d",$5;exit}}' -`

        # Extract the current year of monthly averages
        ccc select ${var}mod datamavgr <<EOF
  SELECT. STEPS${y01}${y12}         1${mlev}${mlev}      ${var}
EOF

        # Create annual average
        timavg datamavgr datamavgr1

        ggstat datamavgr1

        ccc relabl datamavgr1 datamavg  <<EOF
    RELABL GRID               
           GRID                                  0
EOF

        # determine difference from obs data
        sub datamavg ${var}obs diff

        ####### BUG Is this supposed to happen for sss
        if [ 1 -eq 2 -a "$VAR" = "sss" ] ; then
          echo "  XLIN        1000.0       0.0" > IN
          ccc xlin diff diffs input=IN
          echo "  XLIN           1.0       0.0" > IN
          ccc xlin diffs diff input=IN
        fi

        echo "$VAR difference for year $yyy"
        ggstat diff

        # Generate plot frame for ncgm file
        [ "$VAR" = "sst" ] && sstplt $yyy $yyy
        [ "$VAR" = "sss" ] && sssplt $yyy $yyy

      done
    done
  fi

  release TEMPmod
  release SALTmod

  . plot.cdk

end_of_script

. endjcl.cdk

#end_of_job