#!/bin/bash -x
# Shell script to compare observations with simulations using PyNCplot
# 
#   observations: different files from create_OBSnetcdf.py and/or UWyoming_snd_nc.py
#   simulations: multiple outputs from different models and runs
#

# Name of the file with the configuration
if test $1 = '-h'; then
  echo "************************************************"
  echo "***            Script to compare             ***"
  echo "*** sounding and single-station observations ***"
  echo "***         with multiple simulations        ***"
  echo "************************************************"
  echo "obs-sim_Comparison [ConfFile](configuration file)"
else
  rootsh=`pwd`

  configfname=$1

####### ###### ##### #### ### ## #

function uploadvars() {
# Function to upload variables to the system from an ASCII file as:
#   [varname] = [value]
  fileval=$1
  errormsg='ERROR -- error -- ERROR -- error'

  if test ! -f ${fileval}; then
    echo ${errormsg}
    echo "  "${fname}": file '"${fileval}"' does not exist!!"
    exit
  fi

  Nlin=`wc -l ${fileval} | awk '{print $1}'`

  ilin=1
  while test ${ilin} -le ${Nlin}; do
    line=`head -n ${ilin} ${fileval} | tail -n 1`
    varname=`echo ${line} | tr '=' ' ' | awk '{print $1}'`
    value=`echo ${line} | tr '=' ' ' | awk '{print $2}'`
    Lvarname=`expr length ${varname}'0'`

    if test ${Lvarname} -gt 1 && test ! ${varname:0:1} = '#'; then
      export ${varname}=${value}
    fi
    ilin=`expr ${ilin} + 1`
  done  
}

uploadvars ${configfname}
echo "END upload -- end UPLOAD -- END upload -- end UPLOAD"

# files from scratch
if test; then ${fscratch} = 'true'; fscratch=true
else fscratch=false; fi

# figures from scratch
if test; then ${gscratch} = 'true'; gscratch=true
else gscratch=false; fi

function num_hex(){
# Fucntion to transform a number to a hexagonal basis
#   num: number to transform
  fname='num_hex'
  num=$1

  if test ${num} -gt 255; then
    echo ${errormsg}
    echo "  "${fname}": value "${num}" too large !!"
    exit
  fi
  if test ${num} -gt 16; then
    xviA=`expr ${num} / 16`
    numA=`expr ${xviA} '*' 16`
    xviB=`expr ${num} - ${numA}`
  else
    xviA=0
    xviB=${num}
  fi

  if test ${xviA} -gt 9; then
    case ${xviA} in
      10)
        hexA='A'
      ;;
      11)
        hexA='B'
      ;;
      12)
        hexA='C'
      ;;
      13)
        hexA='D'
      ;;
      14)
        hexA='E'
      ;;
      15)
        hexA='F'
      ;;
    esac
  else
    hexA=${xviA}
  fi

  if test ${xviB} -gt 9; then
    case ${xviB} in
      10)
        hexB='A'
      ;;
      11)
        hexB='B'
      ;;
      12)
        hexB='C'
      ;;
      13)
        hexB='D'
      ;;
      14)
        hexB='E'
      ;;
      15)
        hexB='F'
      ;;
    esac
  else
    hexB=${xviB}
  fi

  echo ${hexA}${hexB}
}

function WindRose_plot() {
# Function to plot a Wind rose
#   it: time-step to use from the file
#   sndstid: id of the sounding station
#   sndstn: name of the sounding station to appear in plot ('!' for spaces)
#   timeS: String with the format of the actual time ('!' for spaces)
#   ofigdir: folder for the output of the figure
#   timefS: date of the plot to be used for the figure name
#   filen: file to use
#   figheader: header of the figure's file name

    valuesfig=$1

    it=`echo ${valuesfig} | tr '#' ' ' | awk '{print $1}'`
    sndstid=`echo ${valuesfig} | tr '#' ' ' | awk '{print $2}'`
    sndstn=`echo ${valuesfig} | tr '#' ' ' | awk '{print $3}'`
    timeS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $4}'`
    ofigdir=`echo ${valuesfig} | tr '#' ' ' | awk '{print $5}'`
    timefS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $6}'`
    filen=`echo ${valuesfig} | tr '#' ' ' | awk '{print $7}'`
    figheader=`echo ${valuesfig} | tr '#' ' ' | awk '{print $8}'`

    values='pres|-1;time|'${it}':linepoint;multicol;pres;auto;auto;rainbow;auto:'
    values=${values}${sndstn}'!sounding!WindRose!on!'${timeS}'!local!time:'${kfig}
    values=${values}':cardinals:False:WindRose:True'
    ofign=${ofigdir}'/'${figheader}'_'${sndstid}'_'${timefS}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      echo "  Wind Rose on: "${timeS}
      python ${pyHOME}/drawing.py -o draw_WindRose -S ${values} -v ua,va -f ${filen} >& /dev/null
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_WindRose -S "'"${values}"'" -v ua,va -f ${filen}
        exit
      fi
      convert -trim WindRose.png ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_WindRose -S "'"${values}"'" -v ua,va -f "'"${filen}"'" >> ${ofilefigs}
      echo " " >> ${ofilefigs}
    fi
}

function SkewT_logP_plot() {
# Function to plot a SkewT_logP plot
#   it: time step to use from file
#   sndstid: id of the sounding station
#   sndstn: name of the sounding station to appear in plot ('!' for spaces)
#   timeS: String with the format of the actual time ('!' for spaces)
#   ofigdir: folder for the output of the figure
#   timefS: date of the plot to be used for the figure name
#   filen: file to use
#   figheader: header of the figure's file name

    valuesfig=$1

    it=`echo ${valuesfig} | tr '#' ' ' | awk '{print $1}'`
    sndstid=`echo ${valuesfig} | tr '#' ' ' | awk '{print $2}'`
    sndstn=`echo ${valuesfig} | tr '#' ' ' | awk '{print $3}'`
    timeS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $4}'`
    ofigdir=`echo ${valuesfig} | tr '#' ' ' | awk '{print $5}'`
    timefS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $6}'`
    filen=`echo ${valuesfig} | tr '#' ' ' | awk '{print $7}'`
    figheader=`echo ${valuesfig} | tr '#' ' ' | awk '{print $8}'`
  
    values='time|'${it}',pres|-1:auto:auto:'${sndstn}'!sounding!on!'${timeS}'!local!'
    values=${values}'time:'${kfig}':True'
    ofign=${ofigdir}'/'${figheader}'_'${sndstid}'_'${timefS}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      echo "  Sounding on: "$(echo ${timeS} | tr '!' ' ')
      python ${pyHOME}/drawing.py -o draw_SkewT -S ${values} -v ta,tda,pres -f ${filen} >& /dev/null
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_SkewT -S "'"${values}"'" -v ta,tda,pres -f ${filen}
        exit
      fi
      convert -trim SkewT.png ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_SkewT -S "'"${values}"'" -v ta,tda,pres -f "'"${filen}"'" >> ${ofilefigs}
      echo " " >> ${ofilefigs}
    fi
}

function multi_soundings_plot() {
# Function to plot a multi_sounding figure
#   sndstid: id of the sounding station
#   sndstn: name of the sounding station to appear in plot ('!' for spaces)
#   tatda_evol_labs: labels for the figure
#   tatda_evol_cols: colors of the lines for the figure
#   timeS3: date in the figure
#   ofigdir: folder for the output figure
#   figheader: header of the figure's file name

  valuesfig=$1

  sndstid=`echo ${valuesfig} | tr '@' ' ' | awk '{print $1}'`
  sndstn=`echo ${valuesfig} | tr '@' ' ' | awk '{print $2}'`
  tatda_evol_labs=`echo ${valuesfig} | tr '@' ' ' | awk '{print $3}'`
  tatda_evol_cols=`echo ${valuesfig} | tr '@' ' ' | awk '{print $4}'`
  tatda_evol_files=`echo ${valuesfig} | tr '@' ' ' | awk '{print $5}'`
  timeS3=`echo ${valuesfig} | tr '@' ' ' | awk '{print $6}'`
  ofigdir=`echo ${valuesfig} | tr '@' ' ' | awk '{print $7}'`
  figheader=`echo ${valuesfig} | tr '@' ' ' | awk '{print $8}'`

  tatda_evol_values='auto:auto:multilines!'${tatda_evol_labs}'!'${tatda_evol_cols}'!'
  tatda_evol_values=${tatda_evol_values}'-!,!2:0,auto:'${sndstn}'!sounding!evolution!on!'
  tatda_evol_values=${tatda_evol_values}${timeS3}'!local!time:'${kfig}':True'

  ofign=${ofigdir}'/'${figheader}'_'${stid}'.'${kfig}
  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
  if test ! -f ${ofign}; then
    echo "  Sounding on: "$(echo ${timeS} | tr '!' ' ')
    python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${tatda_evol_values} -f ${tatda_evol_files} >& /dev/null
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S "'"${tatda_evol_values}"'" -f "'"${tatda_evol_files}"'"
      exit
    fi
    convert -trim multi_SkewT.png ${ofign}
    echo "* "${ofign} >> ${ofilefigs}
    echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S "'"${tatda_evol_values}"'" -f "'"${tatda_evol_files}"'" >> ${ofilefigs}
    echo " " >> ${ofilefigs}
  fi
}

function multi_SkewT_logP_plot() {
# Function to plot a multi_SkewT_logP plot with multiple lines
#   sndstid: id of the sounding station
#   sndstn: name of the sounding station to appear in plot ('!' for spaces)
#   complabs: labels for the plot
#   compcols: colors for the plot
#   complines: type of lines for the plot
#   compfiles: files to use for the plot
#   dateS: String with the format of the actual time ('!' for spaces)
#   ofigdir: folder for the output of the figure
#   datefS: date of the plot to be used for the figure name
#   figheader: header of the figure's file name

    valuesfig=$1
    sndstid=`echo ${valuesfig} | tr '@' ' ' | awk '{print $1}'`
    sndstn=`echo ${valuesfig} | tr '@' ' ' | awk '{print $2}'`
    complabs=`echo ${valuesfig} | tr '@' ' ' | awk '{print $3}'`
    compcols=`echo ${valuesfig} | tr '@' ' ' | awk '{print $4}'`
    complines=`echo ${valuesfig} | tr '@' ' ' | awk '{print $5}'`
    compfiles=`echo ${valuesfig} | tr '@' ' ' | awk '{print $6}'`
    dateS=`echo ${valuesfig} | tr '@' ' ' | awk '{print $7}'`
    ofigdir=`echo ${valuesfig} | tr '@' ' ' | awk '{print $8}'`
    datefS=`echo ${valuesfig} | tr '@' ' ' | awk '{print $9}'`
    figheader=`echo ${valuesfig} | tr '@' ' ' | awk '{print $10}'`

    values='auto:auto:multilines!'${complabs}'!'${compcols}'!'
    values=${values}${complines}'!,!2:0,auto:'${ststn}'!sounding!obs!,!sim!'
    values=${values}'comparison!on!'${dateS}'!UTC:png:True'

    ofign=${ofigdir}'/'${figheader}'_'${sndstid}'_'${datefS}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      echo "  Sounding on: "$(echo ${dateS} | tr '!' ' ')
      python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${values} -f ${compfiles} >& /dev/null
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S "'"${values}"'" -f "'"${compfiles}"'"
        exit
      fi
      convert -trim multi_SkewT.png ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S "'"${values}"'" -f "'"${compfiles}"'" >> ${ofilefigs}
      echo " " >> ${ofilefigs}
    fi
}

function Homoeller_SkewT_map_plot() {
# Function to plot a Hovmoeller SkewT_map
#   sndstid: id of the sounding station
#   sndstn: name of the sounding station to appear in plot ('!' for spaces)
#   sndv: variable to plot
#   ofigdir: folder for the output of the figure
#   datefS: date of the plot to be used for the figure name
#   timevals: values to use for the time-axis in the plot
#   simfilen: name of the simulation file at the given point
#   sndfilen: name of the sounding file
#   expl: label of the experiment
#   expn: name of the experiment
#   cbar: color bar
#   figheader: header of the figure's file name

  valuesfig=$1
  sndstid=`echo ${valuesfig} | tr '#' ' ' | awk '{print $1}'`
  sndstn=`echo ${valuesfig} | tr '#' ' ' | awk '{print $2}'`
  sndv=`echo ${valuesfig} | tr '#' ' ' | awk '{print $3}'`
  ofigdir=`echo ${valuesfig} | tr '#' ' ' | awk '{print $4}'`
  datefS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $5}'`
  timevals=`echo ${valuesfig} | tr '#' ' ' | awk '{print $6}'`
  simfilen=`echo ${valuesfig} | tr '#' ' ' | awk '{print $7}'`
  sndfilen=`echo ${valuesfig} | tr '#' ' ' | awk '{print $8}'`
  expl=`echo ${valuesfig} | tr '#' ' ' | awk '{print $9}'`
  expn=`echo ${valuesfig} | tr '#' ' ' | awk '{print $10}'`
  cbar=`echo ${valuesfig} | tr '#' ' ' | awk '{print $11}'`
  figheader=`echo ${valuesfig} | tr '#' ' ' | awk '{print $12}'`

  cfiles=${simfilen}';'${sndv}';time;pres;time|-1,pres|-1@'
  cfiles=${cfiles}${sndfilen}';'${sndv}';time;pres;time|-1,pres|-1'
  cvalues=${sndv}';y;auto|'${timevals}';Vfix,auto,50.,auto'
  cvalues=${cvalues}';'${cbar}',auto,auto;Srange,Srange;auto;obs!,!'${expn}'!'
  cvalues=${cvalues}${sndstn}'!sounding;'${kfig}';flip@y;None;yes'

  ofign=${ofigdir}'/'${figheader}'_'${expl}'_'${sndstid}'_'${sndv}'.'${kfig}
  if ${gscratch}; then rm ${ofign} >& /dev/null; fi
  if test ! -f ${ofign}; then
    python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f ${cfiles} -S ${cvalues}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f "'"${cfiles}"'" -S "'"${cvalues}"'"
      exit
    fi
    convert -trim 2Dshad_obs-sim_comparison_time.png ${ofign}
    echo "* "${ofign} >> ${ofilefigs}
    echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc_time -f "'"${cfiles}"'" -S "'"${cvalues}"'" >> ${ofilefigs}
    echo " " >> ${ofilefigs}
  fi
}

function multiple_time_series_plot(){
# Function to plot multiple time-series
#   sti: id of the surface station
#   stn: name of the surface station to appear in plot ('!' for spaces)
#   sfcv: name of variable to plot
#   ofigdir: folder for the output of the figure
#   timevals: values of the time-axis to use
#   files: files to use
#   labs: labels of the legend
#   lcol: color of lines to use
#   nsfc: minimum value to plot
#   xsfc: maximum value to plot
#   figheader: header of the figure's file name

    valuesfig=$1
    stid=`echo ${valuesfig} | tr '@' ' ' | awk '{print $1}'`
    stn=`echo ${valuesfig} | tr '@' ' ' | awk '{print $2}'`
    sfcv=`echo ${valuesfig} | tr '@' ' ' | awk '{print $3}'`
    ofigdir=`echo ${valuesfig} | tr '@' ' ' | awk '{print $4}'`
    timevals=`echo ${valuesfig} | tr '@' ' ' | awk '{print $5}'`
    files=`echo ${valuesfig} | tr '@' ' ' | awk '{print $6}'`
    labs=`echo ${valuesfig} | tr '@' ' ' | awk '{print $7}'`
    lcol=`echo ${valuesfig} | tr '@' ' ' | awk '{print $8}'`
    nsfc=`echo ${valuesfig} | tr '@' ' ' | awk '{print $9}'`
    xsfc=`echo ${valuesfig} | tr '@' ' ' | awk '{print $10}'`
    figheader=`echo ${valuesfig} | tr '@' ' ' | awk '{print $11}'`

    timen=`echo ${timevals} | tr '|' ' ' | awk '{print $3}'`
    fmtT=`echo ${timevals} | tr '|' ' ' | awk '{print $1"|"$2}'`

    values='ststimes,time;y;'${timen}';auto;'${labs}';'${sfcv}';'${sfcv}
    values=${values}'|evolution|at|'$(echo ${stn} | tr '!' '|')';'${nsfc}','${xsfc}';'
    values=${values}'time|'${fmtT}';0|12;'
    values=${values}${kfig}';-;'${lcol}';.;2.;2.;all;-1;True'

    ofign=${ofigdir}'/'${figheader}'_'${sfcv}_${sti}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv}
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_lines_time -f "'"${files}"'" -S "'"${values}"'" -v ${sfcv}
        exit
      fi
      convert -trim lines_time.${kfig} ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_lines_time -f "'"${files}"'" -S "'"${values}"'" -v ${sfcv} >> ${ofilefigs}
      echo " "${ofign} >> ${ofilefigs}
    fi
}

function shad_contdisc_map_plot(){
# Function to plot a map of a continuos and a discontinuos variable
#   expl: label of the experiment
#   expn: name of the experiment in the title of the plot
#   sfcv: name of variable to plot
#   ofigdir: folder for the output of the figure
#   obsfiles: observtional file to use
#   simfile: simulation file to use
#   it: time-step from simulation
#   oit: time-step from observations
#   timeS: format of time to appear in the title
#   timefS: format of time to be used for the file name
#   nsfc: minimum value to plot
#   xsfc: maximum value to plot
#   mapv: map values
#   mapcover: cover of the map
#   figheader: header of the figure's file name

    valuesfig=$1

    expl=`echo ${valuesfig} | tr '#' ' ' | awk '{print $1}'`
    expn=`echo ${valuesfig} | tr '#' ' ' | awk '{print $2}'`
    sfcv=`echo ${valuesfig} | tr '#' ' ' | awk '{print $3}'`
    ofigdir=`echo ${valuesfig} | tr '#' ' ' | awk '{print $4}'`
    obsfile=`echo ${valuesfig} | tr '#' ' ' | awk '{print $5}'`
    simfile=`echo ${valuesfig} | tr '#' ' ' | awk '{print $6}'`
    it=`echo ${valuesfig} | tr '#' ' ' | awk '{print $7}'`
    oit=`echo ${valuesfig} | tr '#' ' ' | awk '{print $8}'`
    timeS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $9}'`
    timefS=`echo ${valuesfig} | tr '#' ' ' | awk '{print $10}'`
    nsfc=`echo ${valuesfig} | tr '#' ' ' | awk '{print $11}'`
    xsfc=`echo ${valuesfig} | tr '#' ' ' | awk '{print $12}'`
    mapv=`echo ${valuesfig} | tr '#' ' ' | awk '{print $13}'`
    mapcover=`echo ${valuesfig} | tr '#' ' ' | awk '{print $14}'`
    figheader=`echo ${valuesfig} | tr '#' ' ' | awk '{print $15}'`

    CFvarvals=`python $pyHOME/generic.py -o variables_values -S ${sfcv}`
    cbar=`echo ${CFvarvals} | tr ':' ' ' | awk '{print $7}'`

    cfiles=${simfile}';'${sfcv}';XLONG;XLAT;Time|'${it}',time|'${it}','
    cfiles=${cfiles}'west_east|-1,south_north|-1@'${obsfile}';'${sfcv}
    cfiles=${cfiles}';stslon;stslat;time|'${oit}',lon|-1,lat|-1'
    cvalues=${sfcv}':west_east,south_north:auto:'${cbar}',auto,auto:'${nsfc}','
    cvalues=${cvalues}${xsfc}':''auto:obs!,!'${expn}'!'${sfcv}'!on!'${timeS}'!UTC:'
    cvalues=${cvalues}${kfig}':None:'${mapv}':'${mapcover}':yes'

    ofign=${ofigdir}'/'${figheader}'_'${timefS}'_'${sfcv}'_'${expl}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f ${cfiles} -S ${cvalues}
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f "'"${cfiles}"'" -S "'"${cvalues}"'"
      fi
      convert -trim 2Dshad_obs-sim_comparison.png ${ofign} 
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f "'"${cfiles}"'" -S "'"${cvalues}"'" >> ${ofilefigs}
      echo " " >> ${ofilefigs}
    fi
}

#######    #######
## MAIN
    #######

mkdir -p ${odir}

ofilefigs=${wdir}/allins_figures.inf
if ${fscratch}; then 
  echo "*** Instruction from all figures" > ${ofilefigs}
fi

###
## Preparing simulations
###

# sounding diags
diagns=`echo ${wrfsnddiags} | tr ':' ' '`
snddiagd='time@time,bottom_top@P,south_north@XLAT,west_east@XLONG'
snddiagvals=''
idiag=1
for diagn in ${diagns}; do
  diagv=`cat $pyHOME/diagnostics.inf | grep ${diagn}',' | tr ',' ' ' |               \
    awk '{print $2"|"$3}'`
  if test ${idiag} -eq 1; then
    snddiagvals=${diagv}
  else
    snddiagvals=${snddiagvals}','${diagv}
  fi
  idiag=`expr ${idiag} + 1`
done

if test ${sfcvars} != 'None'; then
# surface diags
diagns=`echo ${wrfsfcdiags} | tr ':' ' '`
sfcdiagd='time@time,south_north@XLAT,west_east@XLONG'
sfcdiagvals=''
idiag=1
iadv=1
for diagn in ${diagns}; do
  diagv=`cat $pyHOME/diagnostics.inf | grep WRF | grep ${diagn}',' | tr ',' ' ' |    \
    awk '{print $2"|"$3}'`
  Ldiagv=`expr length $(echo ${diagv} | tr ' ' '@')0`
  if test ${Ldiagv} -lt 2; then
    echo ${errmsg}
    echo "  diagnostic '"${diagn}"' (for WRF) does not exist !!"
    if test ${iadv} -eq 1; then advs=${diagn}; iadv=`expr ${iadv} + 1`
    else advs=${adv}','${diagn}; fi
  else
    if test ${idiag} -eq 1; then
      sfcdiagvals=${diagv}
      idiag=`expr ${idiag} + 1`
    else
      sfcdiagvals=${sfcdiagvals}','${diagv}
    fi
  fi
done
fi

exps=`echo ${sims} | tr ':' ' '`

iexp=1
for exp in ${exps}; do
  expf=`echo ${exp} | tr ',' ' ' | awk '{print $1}'`
  expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
  expn=`echo ${exp} | tr ',' ' ' | awk '{print $3}'`
  expc=`echo ${exp} | tr ',' ' ' | awk '{print $4}'`

  # Joining files
  simjoinselvars=${odir}'/simout_vars_'${expl}'.nc'
  if ${fscratch}; then rm ${simjoinselvars}; fi
  if test ! -f ${simjoinselvars}; then
    values=${simdir}'/'${expf}','${simtimedimn}','${simtimevarn}
    HMT=${simHMT}
    python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values}         \
     -f ${HMT} -v ${simvars}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S "'"${values}"'"  \
        -f "'"${HMT}"'" -v "'"${simvars}"'"
      exit
    fi
    mv netcdf_fold_concatenated_HMT.nc ${simjoinselvars}
  fi

  # Computing sounding diagnostics
  simdiagsf=${odir}'/simout_snddiags_'${expl}'.nc'
  if ${fscratch}; then rm ${simdiagsf}; fi
  if test ! -f ${simdiagsf}; then
    python ${pyHOME}/diagnostics.py -f ${simjoinselvars} -d ${snddiagd} -v ${snddiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simjoinselvars} -d "'"${snddiagd}"'" -v "'"${snddiagvals}"'"
      exit
    fi
    mv diagnostics.nc ${simdiagsf}
    #To deg
    python ${pyHOME}/nc_var.py -o valmod -S subc,273.15 -f ${simdiagsf} -v ta
    #To Pa
    python ${pyHOME}/nc_var.py -o valmod -S divc,100. -f ${simdiagsf} -v pres
  fi

  if test ${sfcvars} != 'None'; then
  # Computing surface diagnostics
  simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'
  if ${fscratch}; then rm ${simsfcdiagf}; fi
  if test ! -f ${simsfcdiagf}; then
    python ${pyHOME}/diagnostics.py -f ${simjoinselvars} -d ${sfcdiagd} -v ${sfcdiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simjoinselvars} -d "'"${sfcdiagd}"'" -v ${sfcdiagvals}
      exit
    fi
    mv diagnostics.nc ${simsfcdiagf}
    #To rad
    python ${pyHOME}/nc_var.py -o valmod -S mulc,0.0174532925199 -f ${simsfcdiagf} -v wds
    python ${pyHOME}/nc_var.py -o varaddattr -S 'units|rad' -f ${simsfcdiagf} -v wds
  fi

  #Adding non-diagnostic variables
  if test ${iadv} -ne 1; then
    varsadd=`echo ${advs} | tr ',' ' '`
    for vadd in ${varsadd}; do
      python ${pyHOME}/nc_var.py -o fvaradd -S ${simjoinselvars},${vadd} -f ${simsfcdiagf}
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simjoinselvars},${vadd}"'" -f ${simsfcdiagf}
        exit
        rm ${simsfcdiagf}
      fi
      # CF varname
      CFvarn=`python $pyHOME/generic.py -o variables_values -S ${vadd} | tr ':' ' ' |  \
        awk '{print $1}'`
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python $pyHOME/generic.py -o variables_values -S ${vadd}
        exit
        rm ${simsfcdiagf}
      fi
      python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${simsfcdiagf} -v ${vadd}
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn} -f ${simsfcdiagf} -v ${vadd}
        exit
        rm ${simsfcdiagf}
      fi
    done
  fi
  fi

  iexp=`expr ${iexp} + 1`

# End experiments
done

mkdir -p ${ofigdir}

###
## Soundings figures
###

# Files in sounding folder (from UWyoming_netcdf.py)
sndfiles=`ls -1 ${snddir}/${sndHf}*`

echo "Soundings files:"${sndfiles}

isnd=1
for sndf in ${sndfiles}; do
  stid=`python $pyHOME/nc_var.py -o grattr -f ${sndf} -S Station_number`

  # Getting wind components
  ouavasndf=${odir}'/sounding_uava_'${stid}'.nc'
  if ${fscratch}; then rm ${ouavasndf}; fi
  if test ! -f ${ouavasndf}; then
    python ${pyHOME}/diagnostics.py -f ${sndf} -d 'time@time,pres@pres' -v 'uavaFROMwswd|ws@wd'
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${sndf} -d 'time@time,pres@pres' -v 'uavaFROMwswd|ws@wd'
      exit
    fi 
    mv diagnostics.nc ${ouavasndf}
  fi
  if test ${isnd} -eq 1; then sndstids=${stid}'@'${sndf}
  else sndstids=${sndstids}':'${stid}'@'${sndf}; fi
  isnd=`expr ${isnd} + 1`
done

sndids=`echo ${sndstids} | tr ':' ' '`

for stidf in ${sndids}; do
  stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
  sndorigfn=`echo ${stidf} | tr '@' ' ' | awk '{print $2}'`
  echo ${stid}" ..."
  it=0
  tatda_evol_values=''
  tatda_evol_files=''
  sndtimes=''
  snddates=''
  sndFdates=''
  snddimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime -f ${sndorigfn} -v time | \
     grep 'dimt:' | awk '{print $2-1}'`

  # Time characteristics
  tunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f ${sndorigfn} | \
    grep -v allvattrs | grep units | awk '{print $3}'`
  Tu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
  Dref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
  Dref=${Dref0:0:4}${Dref0:5:2}${Dref0:8:2}

  # Getting values for evolution plot
  col=`echo ${icolta} | tr ':' ' ' | awk '{print $1}'`
  col2=`echo ${icolta} | tr ':' ' ' | awk '{print $2}'`
  while test ${it} -le ${snddimt}; do
    timev=`python $pyHOME/nc_var.py -o varout -S time:${it} -f ${sndorigfn} -v time | \
      grep NC | awk '{printf ("%d", $2)}'`
    timeS=`date +${figsndTfmt} -d"${Dref} ${timev} ${Tu}"`
    timeS2=`date +%d"$^{"%H"}$" -d"${Dref} ${timev} ${Tu}"`
    timeS3=`date +%Y"-"%b -d"${Dref} ${timev} ${Tu}"`
    timeS4=`date +%Y -d"${Dref} ${timev} ${Tu}"`
    echo "  "${it}" "${timev}" "${Dref}" '"${timev}"' "${Tu}" "$(echo ${timeS} | tr '!' ' ')
    timefS=`date +%Y%m%d%H%M%S -d"${Dref} ${timev} ${Tu}"`

    # for now...
    sndstn=${stid}

    fivals=${it}'#'${stid}'#'${sndstn}'#'${timeS}'#'${ofigdir}'#'${timefS}'#'
    fivals=${fivals}${ouavasndf}'#WindRose'
    WindRose_plot ${fivals}

    fivals=${it}'#'${stid}'#'${sndstn}'#'${timeS}'#'${ofigdir}'#'${timefS}'#'
    fivals=${fivals}${sndorigfn}'#SkewT-logP_obs_ta-tda'
    SkewT_logP_plot ${fivals}

    colH=`num_hex ${col}`
    colH2=`num_hex ${col2}`
    if test ${it} -eq 0; then
      sndtimes=${timev}
      snddates=${timeS}
      sndFdates=${timefS}
      tatda_evol_files=${sndorigfn}':pres|-1,time|'${it}':ta,pres;'
      tatda_evol_files=${tatda_evol_files}${sndorigfn}':pres|-1,time|'${it}':tda,pres'
      tatda_evol_labs='ta_'${timeS2}',tda'
      tatda_evol_cols='#'${colH}${colH2}${colH2}',#'${colH2}${colH2}${colH}
    else
      sndtimes=${sndtimes}','${timev}
      snddates=${snddates}','${timeS}
      sndFdates=${sndFdates}','${timefS}
      tatda_evol_files=${tatda_evol_files}';'${sndorigfn}':pres|-1,time|'${it}':ta,pres'
      tatda_evol_files=${tatda_evol_files}';'${sndorigfn}':pres|-1,time|'${it}':tda,pres'
      if test $(expr ${it} % 2) -eq 0; then 
        tatda_evol_labs=${tatda_evol_labs}','${timeS2}',None'
      else 
        tatda_evol_labs=${tatda_evol_labs}',None,None'
      fi
      tatda_evol_cols=${tatda_evol_cols}',#'${colH}${colH2}${colH2}',#'${colH2}${colH2}${colH}
    fi
    col=`expr ${col} - ${ddcol}`
    col2=`expr ${col2} - ${ddcol}`

    #exit
    it=`expr ${it} + 1`
  done

  fivals=${stid}'@'${sndstn}'@'${tatda_evol_labs}'@'${tatda_evol_cols}'@'
  fivals=${fivals}${tatda_evol_files}'@'${timeS3}'@'${ofigdir}'@SkewT-logP_obs_evol'
  multi_soundings_plot ${fivals}

# End of soundings
done

# snd-sims comparison

for stidf in ${sndids}; do
  stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
  sndorigfn=`echo ${stidf} | tr '@' ' ' | awk '{print $2}'`
  echo "snd: "${stid}

  sndlon=`ncdump -h ${sndorigfn} | grep Station_longitude | awk '{print $3}'`
  sndlat=`ncdump -h ${sndorigfn} | grep Station_latitude | awk '{print $3}'`
  sndts=`echo ${sndtimes} | tr ',' ' '`

  # for now...
  sndstn=${stid}

  # Time characteristics
  tunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f ${sndorigfn} |            \
    grep -v allvattrs | grep units | awk '{print $3}'`
  Tu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
  Dref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
  Dref=${Dref0:0:4}${Dref0:5:2}${Dref0:8:2}

  obsdimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime, -v time -f ${sndorigfn} |  \
    grep dimt | awk '{print $2-1}'`

  # Plotting
  it=0
  while test ${it} -le ${obsdimt}; do

    otimev=`python $pyHOME/nc_var.py -o varout -S time:${it} -f ${sndorigfn} -v time |\
      grep NC | awk '{printf ("%d", $2)}'`
    dateS=`date +${figsndTfmt} -d"${Dref} ${otimev} ${Tu}"`
    timeS2=`date +%d"$^{"%H"}$" -d"${Dref} ${otimev} ${Tu}"`
    timeS3=`date +%Y"-"%b -d"${Dref} ${otimev} ${Tu}"`
    timeS4=`date +%Y -d"${Dref} ${otimev} ${Tu}"`
    datefS=`date +%Y%m%d%H%M%S -d"${Dref} ${otimev} ${Tu}"`

    iexp=1
    for exp in ${exps}; do
      expf=`echo ${exp} | tr ',' ' ' | awk '{print $1}'`
      expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
      expn=`echo ${exp} | tr ',' ' ' | awk '{print $3}'`
      expc=`echo ${exp} | tr ',' ' ' | awk '{print $4}'`

      simjoinselvars=${odir}'/simout_vars_'${expl}'.nc'
      simdiagsf=${odir}'/simout_snddiags_'${expl}'.nc'

      gridsndv=`python ${pyHOME}/nc_var.py -o get_point -f ${simjoinselvars}         \
        -S 'XLONG:XLAT:Time|0' -v ${sndlon}','${sndlat} | tr ' ' '!'`
      xsnd=`echo ${gridsndv} | tr '!' ' ' | awk '{print $1}' | tr ',' ' ' |          \
        awk '{print $1}'`
      ysnd=`echo ${gridsndv} | tr '!' ' ' | awk '{print $1}' | tr ',' ' ' |          \
        awk '{print $2}'`
      dsnd=`echo ${gridsndv} | tr '!' ' ' | awk '{print $2}'`
      if test ${it} -eq 0 && test ${iexp} -eq 1; then
        echo "    Equivalent simulated sounding grid point: "${gridsnd}" distance: "\
          ${dsnd}
      fi

      # Getting values at the sounding point
      simsndptf=${odir}'/simout_vars_sndpt_'${stid}'_'${expl}'.nc'
      if ${fscratch}; then rm ${simsndptf} >& /dev/null; fi
      if test ! -f ${simsndptf}; then
        values='time,0,-1,1@bottom_top,0,-1,1@south_north,'${ysnd}','${ysnd}',1@'
        values=${values}'west_east,'${xsnd}','${xsnd}',1'
        python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${simdiagsf}       \
          -S ${values} -v all
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${simdiagsf}\
            -S "'"${values}"'" -v all
          exit
        fi
        ofilen=`ls -rt1 ${odir} | tail -n 1`
        mv ${odir}/${ofilen} ${simsndptf}
      fi

      # Time information from simulation
      simdimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime, -v time                \
        -f ${simdiagsf} | grep dimt | awk '{print $2}'`
      simtunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f ${simdiagsf} |     \
        grep -v allvattrs | grep units | awk '{print $3}'`
      simTu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
      simDref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
      simDref=${simDref0:0:4}${simDref0:5:2}${simDref0:8:2}

      # Time step from simulated file
      timestepv=`python ${pyHOME}/nc_var.py -o get_time -f ${simsndptf}               \
        -S ${otimev}';'${tunits} -v time | tr ' ' '@'`
      timestep=`echo ${timestepv} | tr '@' ' ' | awk '{print $1}'`
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/nc_var.py -o get_time -f ${simsndptf}                  \
          -S ${otimev}';'${tunits} -v time
        exit
      fi

      if test ${iexp} -eq 1; then
        echo "     "${it}": it="${timestep}" "$(echo ${dateS} | tr '!' ' ')" "${datefS}
        tacomplabs='$ta^{obs}$,$ta^{'${expn}'}$'
        tdacomplabs='$tda^{obs}$,$tda^{'${expn}'}$'
        tacompcols='#FF6464,'${expc}
        tdacompcols='#6464FF,'${expc}
        tacompltyp='-,-'
        tdacompltyp='-,-.'
        tacompfiles=${sndorigfn}':pres|-1,time|'${it}':ta,pres;'
        tacompfiles=${tacompfiles}${simsndptf}':bottom_top|-1,time|'${timestep}':'
        tacompfiles=${tacompfiles}'ta,pres'
        tdacompfiles=${sndorigfn}':pres|-1,time|'${it}':tda,pres;'
        tdacompfiles=${tdacompfiles}${simsndptf}':bottom_top|-1,time|'${timestep}':'
        tdacompfiles=${tdacompfiles}'tda,pres'
      else
        tacomplabs=${tacomplabs}',$ta^{'${expn}'}$'
        tdacomplabs=${tdacomplabs}',$tda^{'${expn}'}$'
        tacompcols=${tacompcols}','${expc}
        tdacompcols=${tdacompcols}','${expc}
        tacompltyp=${tacompltyp}',-'
        tdacompltyp=${tdacompltyp}',-.'
        tacompfiles=${tacompfiles}';'${simsndptf}':bottom_top|-1,time|'${timestep}':'
        tacompfiles=${tacompfiles}'ta,pres'
        tdacompfiles=${tdacompfiles}';'${simsndptf}':bottom_top|-1,time|'${timestep}
        tdacompfiles=${tdacompfiles}':tda,pres'
      fi

      # Hovmoeller maps
      sndvs=`echo ${sndvars} | tr ':' ' '`
      for sndv in ${sndvs}; do

        CFvarvals=`python $pyHOME/generic.py -o variables_values -S ${sndv}`
        cbar=`echo ${CFvarvals} | tr ':' ' ' | awk '{print $7}'`

        fivals=${stid}'#'${sndstn}'#'${sndv}'#'${ofigdir}'#'${datefS}'#'${fmtTts}
        fivals=${fivals}'#'${simsndptf}'#'${sndorigfn}'#'${expl}'#'${expn}'#'${cbar}
        fivals=${fivals}'#SkewT-logP_obs-sim_evol'

        Homoeller_SkewT_map_plot ${fivals}

      # end sounding vars
      done

      iexp=`expr ${iexp} + 1`
    # End experiments
    done

    complabs=${tacomplabs}','${tdacomplabs}
    compcols=${tacompcols}','${tdacompcols}
    compltyp=${tacompltyp}','${tdacompltyp}
    compfiles=${tacompfiles}';'${tdacompfiles}

    fivals=${stid}'@'${sndstn}'@'${complabs}'@'${compcols}'@'${compltyp}'@'
    fivals=${fivals}${compfiles}'@'${dateS}'@'${ofigdir}'@'${datefS}'@'${ofigdir}'@'
    fivals=${fivals}'SkewT-logP_obs-sim_step'

    multi_SkewT_logP_plot ${fivals}

    # Wind Rose comparison
    fivals=${it}'#'${stid}'#'${sndstn}'#'${timeS}'#'${ofigdir}'#'${timefS}'#'
    fivals=${fivals}${ouavasndf}'#WindRose_obs-sim_'
    multi_WindRose_plot ${fivals}

    it=`expr ${it} + 1`
  # End sounding time-steps
  done

# End of soundings stations
done

###
## SFC figures
###

# Joining all observations
obssfcfile=${odir}'/all_single-stations.nc'
if test ! -f ${wdir}/${sfcobsdir}/${obssfcfile}; then
  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles -S                       \
    ${sfcobsdir}':OBSnetcdf' -v all
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles \
      -S ${sfcobsdir}':OBSnetcdf' -v all
  fi
  mv joined_singlestations.nc ${obssfcfile}
fi

# getting min,max values
if test ${sfcvars} != 'None'; then
fcs=`echo ${sfcvars} | tr ':' ' '`

iv=1
for sfcpv in ${sfcs}; do
  sfcv=`echo ${sfcpv} | tr '@' ' ' | awk '{print $1}'`
  nsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $2}'`
  xsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $3}'`
  if test ${nsfc} = 'FROMobs'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc=`echo ${fstats} | awk '{print $3}'`
  elif test ${nsfc} = 'FROMsims'; then
    iexp=1
    for exp in ${exps}; do
      expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
      simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'
    
      fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None            \
        -f ${simsfcdiagf} -v ${sfcv} | grep MAT | grep ${sfcv}`
      if test ${iexp} -eq 1; then
        nsfc=`echo ${fstats} | awk '{print $3}'`
      else
        nsfc2=`echo ${fstats} | awk '{print $3}'`
        nsfcf=`echo ${nsfc}" "${nsfc2} | awk '{if ($1 < $2) {print $1;} else {print $2;}}'`
      fi
    done
  elif test ${nsfc} = 'FROMobssims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc=`echo ${fstats} | awk '{print $3}'`
    for exp in ${exps}; do
      expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
      simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'
    
      fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None            \
        -f ${simsfcdiagf} -v ${sfcv} | grep MAT | grep ${sfcv}`
      nsfc2=`echo ${fstats} | awk '{print $3}'`
      nsfc=`echo ${nsfc}" "${nsfc2} | awk '{if ($1 < $2) {print $1;} else {print $2;}}'`
    done
  fi
  if test ${xsfc} = 'FROMobs'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    xsfc=`echo ${fstats} | awk '{print $4}'`
  elif test ${xsfc} = 'FROMsims'; then
    iexp=1
    for exp in ${exps}; do
      expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
      simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'
    
      fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None            \
        -f ${simsfcdiagf} -v ${sfcv} | grep MAT | grep ${sfcv}`
      if test ${iexp} -eq 1; then
        xsfc=`echo ${fstats} | awk '{print $4}'`
      else
        xsfc2=`echo ${fstats} | awk '{print $4}'`
        xsfcf=`echo ${xsfc}" "${xsfc2} | awk '{if ($1 > $2) {print $1;} else {print $2;}}'`
      fi
    done
  elif test ${nsfc} = 'FROMobssims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    xsfc=`echo ${fstats} | awk '{print $4}'`
    iexp=1
    for exp in ${exps}; do
      expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
      simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'
    
      fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None            \
        -f ${simsfcdiagf} -v ${sfcv} | grep MAT | grep ${sfcv}`
      xsfc2=`echo ${fstats} | awk '{print $4}'`
      xsfcf=`echo ${xsfc}" "${xsfc2} | awk '{if ($1 > $2) {print $1;} else {print $2;}}'`
    done
  fi
  if test $iv -eq 1; then
    newsfcvars=${sfcv}'@'${nsfc}'@'${xsfc}
  else
    newsfcvars=${newsfcvars}':'${sfcv}'@'${nsfc}'@'${xsfc}
  fi
  iv=`expr ${iv} + 1`
done
sfcs=`echo ${newsfcvars} | tr ':' ' '`
echo "new sfc variables: "${sfcs}

# Time-series

# stations
# FROM: https://stackoverflow.com/questions/2961635/using-awk-to-print-all-columns-from-the-nth-to-the-last
stsids=`python ${pyHOME}/nc_var.py -o varout -v id -f \
  ${obssfcfile} -S 'Nstations:-1|time:0' | awk '{$1=""; printf("%d:",substr($0,2))}'`
stsnames=`python ${pyHOME}/nc_var.py -o varout -v stsname -f \
  ${obssfcfile} -S Nstations:-1 | awk '{$1=""; print substr($0,2)}' | tr '\n' ':' | \
  tr ' ' '!'`
stslons=`python ${pyHOME}/nc_var.py -o varout -v stslon -f \
  ${obssfcfile} -S Nstations:-1 | awk '{$1=""; print substr($0,2)}' | tr '\n' ':'`
stslats=`python ${pyHOME}/nc_var.py -o varout -v stslat -f \
  ${obssfcfile} -S Nstations:-1 | awk '{$1=""; print substr($0,2)}' | tr '\n' ':'`

Nsts=`echo ${stsnames} | tr ':' ' ' | wc -w | awk '{print $1}'`
ist=1
while test ${ist} -le ${Nsts}; do
  sti=`echo ${stsids} | tr ':' '\n' | head -n ${ist} | tail -n 1`
  stn=`echo ${stsnames} | tr ':' '\n' | head -n ${ist} | tail -n 1`
  stl=`echo ${stslons} | tr ':' '\n' | head -n ${ist} | tail -n 1`
  stL=`echo ${stslats} | tr ':' '\n' | head -n ${ist} | tail -n 1`
  ist_1=`expr ${ist} - 1`

  sfccompfiles=${obssfcfile}'%time|-1;Nstations|'${ist_1}
  sfccomplabs='obs'
  sfccompcols='#000000'

  iexp=1
  for exp in ${exps}; do
    expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
    expn=`echo ${exp} | tr ',' ' ' | awk '{print $3}'`
    expc=`echo ${exp} | tr ',' ' ' | awk '{print $4}'`

    simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'

    gridsfcv=`python ${pyHOME}/nc_var.py -o get_point -f ${simsfcdiagf}              \
      -S 'XLONG:XLAT:Time|0' -v ${stl}','${stL} | tr ' ' '!'`
    xgrid=`echo ${gridsfcv} | tr '!' ' ' | awk '{print $1}' | tr ',' ' ' |           \
      awk '{print $1}'`
    ygrid=`echo ${gridsfcv} | tr '!' ' ' | awk '{print $1}' | tr ',' ' ' |           \
      awk '{print $2}'`
    distgrid=`echo ${gridsfcv} | tr '!' ' ' | awk '{print $2}'`

    sfcstats=${sti}'@'${stn}'@'${stl}'@'${stL}'@'${gridsfc}

    sfccompfiles=${sfccompfiles}','${simsfcdiagf}'%time|-1;south_north|'${ygrid}';west_east|'
    sfccompfiles=${sfccompfiles}${xgrid}
    sfccomplabs=${sfccomplabs}','${expn}
    sfccompcols=${sfccompcols}','${expc}
  done
  # end experiments

  for sfcpv in ${sfcs}; do
    sfcv=`echo ${sfcpv} | tr '@' ' ' | awk '{print $1}'`
    nsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $2}'`
    xsfc=`echo ${sfcpv} | tr '@' ' ' | awk '{print $3}'`

    fivals=${sti}'@'${stn}'@'${sfcv}'@'${ofigdir}'@'${fmtTts}'@'${sfccompfiles}'@'
    fivals=${fivals}${sfccomplabs}'@'${sfccompcols}'@'${nsfc}'@'${xsfc}
    fivals=${fivals}'@obs-sim_evol_ts'

    multiple_time_series_plot ${fivals}

  done
  # end observations

  ist=`expr $ist + 1`
done
# end sfc statinos

# Cont-disc comparison
iexp=1
for exp in ${exps}; do
  expl=`echo ${exp} | tr ',' ' ' | awk '{print $2}'`
  expn=`echo ${exp} | tr ',' ' ' | awk '{print $3}'`
  expc=`echo ${exp} | tr ',' ' ' | awk '{print $4}'`

  simsfcdiagf=${odir}'/simout_sfcdiags_'${expl}'.nc'

  # Time information from simulation
  simdimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime, -v time -f ${simsfcdiagf}  \
    | grep dimt | awk '{print $2-1}'`
  simtunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f ${simsfcdiagf} |       \
    grep -v allvattrs | grep units | awk '{print $3}'`
  simTu=`echo ${tunits} | tr '!' ' ' | awk '{print $1}'`
  simDref0=`echo ${tunits} | tr '!' ' ' | awk '{print $3}'`
  simDref=${simDref0:0:4}${simDref0:5:2}${simDref0:8:2}

  it=0
  while test ${it} -le ${simdimt}; do
    simtstep=`python ${pyHOME}/nc_var.py -o varout -f ${simsfcdiagf} -S              \
      'time:'${it} -v time | grep NC | awk '{printf ("%d",$2)}'`
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o varout -f ${simsfcdiagf} -S 'time:'${it} -v time
      exit
    fi

    timeS=`date +${figsndTfmt} -d"${simDref} ${simtstep} ${simTu}"`
    echo ${it}" "${simtstep}" "${simDref}" "${simTu}" "$(echo ${timeS} | tr '!' ' ')
    timefS=`date +%Y%m%d%H%M%S -d"${simDref} ${simtstep} ${simTu}"`

    # Looking for obs time
    oitv=`python ${pyHOME}/nc_var.py -o get_time -f ${obssfcfile}                    \
      -S ${simtstep}';'${tunits} -v ststimes | grep -v differ | tr ' ' '@'`
    oit=`echo ${oitv} | tr '@' ' ' | awk '{print $1}'`

    for sfcnxv in ${sfcs}; do
      sfcv=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $1}'`
      nsfc=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $2}'`
      xsfc=`echo ${sfcnxv} | tr '@' ' ' | awk '{print $3}'`

      fivals=${expl}'#'${expn}'#'${sfcv}'#'${ofigdir}'#'${obssfcfile}'#'${simsfcdiagf}
      fivals=${fivals}'#'${it}'#'${oit}'#'${timeS}'#'${timefS}'#'${nsfc}'#'${xsfc}'#'
      fivals=${fivals}${mapv}'#'${mapcover}'#map_obs-sim'

      shad_contdisc_map_plot ${fivals}

    # end of variables
    done

  it=`expr ${it} + 1`
  # end of time-steps
  done

# end of simulations
done
fi

fi
