#!/bin/bash
# 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 ${fscratch} = 'true'; then fscratch=true
else fscratch=false; fi

# figures from scratch
if test ${gscratch} = 'true'; then 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 isin_list() {
# Function to check whether a value is in a list
  list=$1
  value=$2
  
  is=`echo ${list} | tr ':' '\n' | awk '{print "@"$1"@"}' | grep '@'${value}'@' | wc -w`
  if test ${is} -eq 1
  then
    true
  else
    false
  fi
}

function get_WRFsounding_data(){
# Function to get the data to compute analysis in a sounding point
#   stlon= longitude value of the station
#   stlat= latitude value of the station
#   filen= name of the file
#   lonvarn= name of the longitude variable
#   latvarn= name of the latitude variable
#   snddiagnames= diagnostics to compute (as ':' separated list)
#   snddiagd= name of the dimensions and dimension variables to use for diagnostics

  funcvals=$1

  stlon=`echo ${funcvals} | tr '#' ' ' | awk '{print $1}'`
  stlat=`echo ${funcvals} | tr '#' ' ' | awk '{print $2}'`
  filen=`echo ${funcvals} | tr '#' ' ' | awk '{print $3}'`
  lonvarn=`echo ${funcvals} | tr '#' ' ' | awk '{print $4}'`
  latvarn=`echo ${funcvals} | tr '#' ' ' | awk '{print $5}'`
  snddiagnames=`echo ${funcvals} | tr '#' ' ' | awk '{print $6}'`
  snddiagd=`echo ${funcvals} | tr '#' ' ' | awk '{print $7}'`

  vals=${lonvarn}':'${latvarn}':Time|0'
  gridsndv=`python ${pyHOME}/nc_var.py -o get_point -f ${filen} -S ${vals}           \
    -v ${stlon},${stlat} | tr ' ' '!'`
  echo "* get grid point: "${stlon}", "${stlat} >> ${ofileins}
  echo python ${pyHOME}/nc_var.py -o get_point -f ${filen} -S "'"${vals}"'"          \
    -v ${stlon},${stlat} >> ${ofileins}
  echo " " >> ${ofileins}

  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}'`

  xsnd1=`expr ${xsnd} + 1`
  ysnd1=`expr ${ysnd} + 1`

  # Getting values
  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@south_north_stag,'${ysnd}','
    values=${values}${ysnd1}',1@west_east_stag,'${xsnd}','${xsnd1}',1:False'
    python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}               \
      -S ${values} -v all
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}        \
        -S "'"${values}"'" -v all
      exit
    fi
    echo "* get values at snd point: "${stlon}", "${stlat} >> ${ofileins}
    echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}          \
      -S "'"${values}"'" -v all >> ${ofileins}
    echo " " >> ${ofileins}

    yslc='south_north_B'${ysnd}'-E'${ysnd}'-I1'
    ystslc='south_north_stag_B'${ysnd}'-E'${ysnd1}'-I1'
    xslc='west_east_B'${xsnd}'-E'${xsnd}'-I1'
    xstslc='west_east_stag_B'${xsnd}'-E'${xsnd1}'-I1'

    Hfilen=`basename ${filen}`
    Dfilen=`dirname ${filen}`
    Hhfilen=`echo ${Hfilen} | tr '.' ' ' | awk '{print $1}'`

    ofilen='newfile_multidims.nc'
    mv ${ofilen} ${simsndptf}
  fi
    
  # Computing sounding diagnostics
  diagns=`echo ${snddiagnames} | tr ':' ' '`

  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 ${fscratch}; then rm ${simsnddiagsf}; fi
  if test ! -f ${simsnddiagsf}; then
    python ${pyHOME}/diagnostics.py -f ${simsndptf} -d ${snddiagd} -v ${snddiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simsndptf} -d "'"${snddiagd}"'" -v   \
        "'"${snddiagvals}"'"
      exit
    fi
    echo "* computing diagnostics at snd point: "${stlon}", "${stlat} >> ${ofileins}
    echo python ${pyHOME}/diagnostics.py -f "'"${simsndptf}"'" -d "'"${snddiagd}"'"  \
      -v "'"${snddiagvals}"'" >> ${ofileins}
    echo " " >> ${ofileins}

    mv diagnostics.nc ${simsnddiagsf}
    if $(isin_list ${snddiagnames} WRFt); then
      #To deg
      python ${pyHOME}/nc_var.py -o valmod -S subc,273.15 -f ${simsnddiagsf} -v ta
    fi
    if $(isin_list ${snddiagnames} WRFp); then
      #To Pa
      python ${pyHOME}/nc_var.py -o valmod -S divc,100. -f ${simsnddiagsf} -v pres
    fi
  fi
}

function get_WRFsurface_data(){
# Function to get the data to compute analysis in a surface point
#   stlon= longitude value of the station
#   stlat= latitude value of the station
#   filen= name of the file
#   lonvarn= name of the longitude variable
#   latvarn= name of the latitude variable
#   sfcdiagnames= diagnostics to compute (as ':' separated list)
#   nondiagvs= ':' list of variables which are not from a diagnostics ('None' for 
#     without)
#   sfcdiagd= name of the dimensions and dimension variables to use for diagnostics

  funcvals=$1

  stlon=`echo ${funcvals} | tr '#' ' ' | awk '{print $1}'`
  stlat=`echo ${funcvals} | tr '#' ' ' | awk '{print $2}'`
  filen=`echo ${funcvals} | tr '#' ' ' | awk '{print $3}'`
  lonvarn=`echo ${funcvals} | tr '#' ' ' | awk '{print $4}'`
  latvarn=`echo ${funcvals} | tr '#' ' ' | awk '{print $5}'`
  sfcdiagnames=`echo ${funcvals} | tr '#' ' ' | awk '{print $6}'`
  nondiagvs=`echo ${funcvals} | tr '#' ' ' | awk '{print $7}'`
  sfcdiagd=`echo ${funcvals} | tr '#' ' ' | awk '{print $8}'`

  vals=${lonvarn}':'${latvarn}':Time|0'
  gridsndv=`python ${pyHOME}/nc_var.py -o get_point -f ${filen} -S ${vals}           \
    -v ${stlon},${stlat} | tr ' ' '!'`
  echo "* get grid point: "${stlon}", "${stlat} >> ${ofileins}
  echo python ${pyHOME}/nc_var.py -o get_point -f ${filen} -S "'"${vals}"'"          \
    -v ${stlon},${stlat} >> ${ofileins}
  echo " " >> ${ofileins}
  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}'`

  xsnd1=`expr ${xsnd} + 1`
  ysnd1=`expr ${ysnd} + 1`

  # Getting values
  if ${fscratch}; then rm ${simsfcptf} >& /dev/null; fi
  if test ! -f ${simsfcptf}; then
    values='Time,0,-1,1@south_north,'${ysnd}','${ysnd}',1@'
    values=${values}'west_east,'${xsnd}','${xsnd}',1@south_north_stag,'${ysnd}','
    values=${values}${ysnd1}',1@west_east_stag,'${xsnd}','${xsnd1}',1:False'
    python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}               \
      -S ${values} -v all
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}        \
        -S "'"${values}"'" -v all
      exit
    fi
  echo "* get values at sfc point: "${stlon}", "${stlat} >> ${ofileins}
  echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}            \
      -S "'"${values}"'" -v all >> ${ofileins}
  echo " " >> ${ofileins}
    yslc='south_north_B'${ysnd}'-E'${ysnd}'-I1'
    ystslc='south_north_stag_B'${ysnd}'-E'${ysnd1}'-I1'
    xslc='west_east_B'${xsnd}'-E'${xsnd}'-I1'
    xstslc='west_east_stag_B'${xsnd}'-E'${xsnd1}'-I1'

    ofilen='newfile_multidims.nc'
    mv ${ofilen} ${simsfcptf}
  fi
    
  # Computing sfc diagnostics
  diagns=`echo ${sfcdiagnames} | tr ':' ' '`

  idiag=1
  for diagn in ${diagns}; do
    diagv=`cat $pyHOME/diagnostics.inf | grep WRF | grep ${diagn}',' | tr ',' ' ' |  \
      awk '{print $2"|"$3}'`
    if test ${idiag} -eq 1; then
      sfcdiagvals=${diagv}
    else
      sfcdiagvals=${sfcdiagvals}','${diagv}
    fi
    idiag=`expr ${idiag} + 1`
  done

  if ${fscratch}; then rm ${simsfcdiagsf}; fi
  if test ! -f ${simsfcdiagsf}; then
    python ${pyHOME}/diagnostics.py -f ${simsfcptf} -d ${sfcdiagd} -v ${sfcdiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simsfcptf} -d "'"${sfcdiagd}"'"      \
        -v "'"${sfcdiagvals}"'"
      exit
    fi
    echo "* diagnostics at sfc point: "${stlon}", "${stlat} >> ${ofileins}
    echo python ${pyHOME}/diagnostics.py -f ${simsfcptf} -d "'"${sfcdiagd}"'"        \
      -v "'"${sfcdiagvals}"'" >> ${ofileins}
    echo " " >> ${ofileins}
    mv diagnostics.nc ${simsfcdiagsf}
    if $(isin_list ${sfcdiagnames} wds); then
      #To rad
      python ${pyHOME}/nc_var.py -o valmod -S mulc,0.0174532925199 -f ${simsfcdiagsf}\
        -v wds
      python ${pyHOME}/nc_var.py -o varaddattr -S 'units|rad' -f ${simsfcdiagsf} -v wds
    fi

    #Adding non-diagnostic variables
    if test ! ${nondiagvs} = 'None'; then
      varsadd=`echo ${nondiagvs} | tr ':' ' '`
      for vadd in ${varsadd}; do
        python ${pyHOME}/nc_var.py -o fvaradd -S ${simsfcptf},${vadd}                \
          -f ${simsfcdiagsf}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcptf},${vadd}"'"   \
            -f ${simsfcdiagsf}
          exit
          rm ${simsfcdiagsf}
        fi
        echo "* adding variable '"${vadd}"' at sfc point: "${stlon}", "${stlat} >>   \
          ${ofileins}
        echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcptf},${vadd}"'"     \
          -f ${simsfcdiagsf} >> ${ofileins}
        echo " " >> ${ofileins}
        # 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 ${simsfcdiagsf}
        fi
        python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${simsfcdiagsf}        \
          -v ${vadd}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn}                  \
            -f ${simsfcdiagsf} -v ${vadd}
          exit
          rm ${simsfcdiagsf}
        fi
      done
    fi
  fi
}

function get_WRFsfc2D_data(){
# Function to get the data to compute analysis in a surface map
#   filen= name of the file
#   sfcdiagnames= diagnostics to compute (as ':' separated list)
#   nondiagvs= ':' list of variables which are not from a diagnostics ('None' for 
#     without)
#   sfcmapdiagd= name of the dimensions and dimension variables to use for diagnostics

  funcvals=$1

  filen=`echo ${funcvals} | tr '#' ' ' | awk '{print $1}'`
  sfcdiagnames=`echo ${funcvals} | tr '#' ' ' | awk '{print $2}'`
  nondiagvs=`echo ${funcvals} | tr '#' ' ' | awk '{print $3}'`
  sfcmapdiagd=`echo ${funcvals} | tr '#' ' ' | awk '{print $4}'`

  # Computing sfc diagnostics
  diagns=`echo ${sfcdiagnames} | tr ':' ' '`

  # getting that variables only required for the diagnostics
  vardiags=''
  idiag=1
  for diagn in ${diagns}; do
    diagv=`cat $pyHOME/diagnostics.inf | grep WRF | grep ${diagn}',' | tr ',' ' ' |  \
      awk '{print $2"|"$3}'`
    if test ${idiag} -eq 1; then
      sfcdiagvals=${diagv}
      vardiags=`echo ${diagv} | tr '|' ' ' | awk '{print $2}' | tr '@' ':'`
    else
      sfcdiagvals=${sfcdiagvals}','${diagv}
      vardiag=`echo ${diagv} | tr '|' ' ' | awk '{print $2}' | tr '@' ' '`
      for vd in ${vardiag}; do
        if ! $(isin_list ${vardiags} ${vd}); then vardiags=${vardiags}':'${vd}; fi
      done
    fi
    idiag=`expr ${idiag} + 1`
  done

  # Inlcuding non-diagnostics
  if test ! ${nondiagvs} = 'None'; then
    nvns=`echo ${nondiagvs} | tr ':' ' '`
    for vn in ${nvns}; do
      if ! $(isin_list ${vardiags} ${vn}); then vardiags=${vardiags}':'${vn}; fi
    done
    if ! $(isin_list ${vardiags} 'time'); then vardiags=${vardiags}':time'; fi
  fi
  vardiag=`echo ${vardiags} | tr ':' ','`
  
  # Getting values
  if ${fscratch}; then rm ${simsfcmapf} >& /dev/null; fi
  if test ! -f ${simsfcmapf}; then
    values='Time,0,-1,1@south_north,0,-1,1@west_east,0,-1,1@south_north_stag,0,-1,1@'
    values=${values}'west_east_stag,0,-1,1:False'
    python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}               \
      -S ${values} -v ${vardiag}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}        \
        -S "'"${values}"'"  -v ${vardiag}
      exit
    fi
    echo "* get values for 2D map diagnostics: " >> ${ofileins}
    echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}          \
      -S "'"${values}"'" -v ${vardiag} >> ${ofileins}
    echo " " >> ${ofileins}

    ofilen='newfile_multidims.nc'
    mv ${ofilen} ${simsfcmapf}
  fi

  if ${fscratch}; then rm ${simsfcmapdiagsf}; fi
  if test ! -f ${simsfcmapdiagsf}; then
    python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d ${sfcmapdiagd}               \
      -v ${sfcdiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d "'"${sfcmapdiagd}"'"  \
        -v "'"${sfcdiagvals}"'"
      exit
    fi
    echo "* diagnostics for sfc 2D map: " >> ${ofileins}
    echo python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d "'"${sfcmapdiagd}"'"    \
      -v "'"${sfcdiagvals}"'" >> ${ofileins}
    echo " " >> ${ofileins}
    mv diagnostics.nc ${simsfcmapdiagsf}
    if $(isin_list ${sfcdiagnames} wds); then
      #To rad
      python ${pyHOME}/nc_var.py -o valmod -S mulc,0.0174532925199                   \
        -f ${simsfcmapdiagsf} -v wds
      python ${pyHOME}/nc_var.py -o varaddattr -S 'units|rad' -f ${simsfcmapdiagsf}  \
        -v wds
    fi

    #Adding non-diagnostic variables
    if test ! ${nondiagvs} = 'None'; then
      varsadd=`echo ${nondiagvs} | tr ':' ' '`
      for vadd in ${varsadd}; do
        python ${pyHOME}/nc_var.py -o fvaradd -S ${simsfcmapf},${vadd}               \
          -f ${simsfcmapdiagsf}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcmapf},${vadd}"'"  \
            -f ${simsfcmapdiagsf}
          exit
          rm ${simsfcmapdiagsf}
        fi
        echo "* adding variable '"${vadd}"' at sfc map: " >> ${ofileins}
        echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcmapf},${vadd}"'"    \
          -f ${simsfcmapdiagsf} >> ${ofileins}
        echo " " >> ${ofileins}
        # 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 ${simsfcmapdiagsf}
        fi
        python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${simsfcmapdiagsf}     \
          -v ${vadd}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn}                  \
            -f ${simsfcmapdiagsf} -v ${vadd}
          exit
          rm ${simsfcmapdiagsf}
        fi
      done
    fi
  fi
}

function get_WRFCDXsfc2D_data(){
# Function to get the data from WRFCDX to compute analysis in a surface map
#   filen= name of the file
#   sfcdiagnames= diagnostics to compute (as ':' separated list)
#   nondiagvs= ':' list of variables which are not from a diagnostics ('None' for 
#     without)
#   sfcmapdiagd= name of the dimensions and dimension variables to use for diagnostics

  funcvals=$1

  filen=`echo ${funcvals} | tr '#' ' ' | awk '{print $1}'`
  sfcdiagnames=`echo ${funcvals} | tr '#' ' ' | awk '{print $2}'`
  nondiagvs=`echo ${funcvals} | tr '#' ' ' | awk '{print $3}'`
  sfcmapdiagd=`echo ${funcvals} | tr '#' ' ' | awk '{print $4}'`

  # Computing sfc diagnostics
  diagns=`echo ${sfcdiagnames} | tr ':' ' '`

  # getting that variables only required for the diagnostics
  vardiags=''
  idiag=1
  for diagn in ${diagns}; do
    diagv=`cat $pyHOME/diagnostics.inf | grep WRF | grep ${diagn}',' | tr ',' ' ' |  \
      awk '{print $2"|"$3}'`
    if test ${idiag} -eq 1; then
      sfcdiagvals=${diagv}
      vardiags=`echo ${diagv} | tr '|' ' ' | awk '{print $2}' | tr '@' ':'`
    else
      sfcdiagvals=${sfcdiagvals}','${diagv}
      vardiag=`echo ${diagv} | tr '|' ' ' | awk '{print $2}' | tr '@' ' '`
      for vd in ${vardiag}; do
        if ! $(isin_list ${vardiags} ${vd}); then vardiags=${vardiags}':'${vd}; fi
      done
    fi
    idiag=`expr ${idiag} + 1`
  done

  # Inlcuding non-diagnostics
  if test ! ${nondiagvs} = 'None'; then
    nvns=`echo ${nondiagvs} | tr ':' ' '`
    for vn in ${nvns}; do
      if ! $(isin_list ${vardiags} ${vn}); then vardiags=${vardiags}':'${vn}; fi
    done
    if ! $(isin_list ${vardiags} 'time'); then vardiags=${vardiags}':time'; fi
  fi
  vardiag=`echo ${vardiags} | tr ':' ','`
  
  # Getting values
  if ${fscratch}; then rm ${simsfcmapf} >& /dev/null; fi
  if test ! -f ${simsfcmapf}; then
    values='Time,0,-1,1@south_north,0,-1,1@west_east,0,-1,1:False'
    python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}               \
      -S ${values} -v ${vardiag}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}        \
        -S "'"${values}"'"  -v ${vardiag}
      exit
    fi
    echo "* get values for 2D map diagnostics: " >> ${ofileins}
    echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${filen}          \
      -S "'"${values}"'" -v ${vardiag} >> ${ofileins}
    echo " " >> ${ofileins}

    ofilen='newfile_multidims.nc'
    mv ${ofilen} ${simsfcmapf}
  fi

  if ${fscratch}; then rm ${simsfcmapdiagsf}; fi
  if test ! -f ${simsfcmapdiagsf}; then
    python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d ${sfcmapdiagd}               \
      -v ${sfcdiagvals}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d "'"${sfcmapdiagd}"'"  \
        -v "'"${sfcdiagvals}"'"
      exit
    fi
    echo "* diagnostics for sfc 2D map: " >> ${ofileins}
    echo python ${pyHOME}/diagnostics.py -f ${simsfcmapf} -d "'"${sfcmapdiagd}"'"    \
      -v "'"${sfcdiagvals}"'" >> ${ofileins}
    echo " " >> ${ofileins}
    mv diagnostics.nc ${simsfcmapdiagsf}

    #Adding non-diagnostic variables
    if test ! ${nondiagvs} = 'None'; then
      varsadd=`echo ${nondiagvs} | tr ':' ' '`
      for vadd in ${varsadd}; do
        python ${pyHOME}/nc_var.py -o fvaradd -S ${simsfcmapf},${vadd}               \
          -f ${simsfcmapdiagsf}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcmapf},${vadd}"'"  \
            -f ${simsfcmapdiagsf}
          exit
          rm ${simsfcmapdiagsf}
        fi
        echo "* adding variable '"${vadd}"' at sfc map: " >> ${ofileins}
        echo python ${pyHOME}/nc_var.py -o fvaradd -S "'"${simsfcmapf},${vadd}"'"    \
          -f ${simsfcmapdiagsf} >> ${ofileins}
        echo " " >> ${ofileins}
        # 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 ${simsfcmapdiagsf}
        fi
        python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${simsfcmapdiagsf}     \
          -v ${vadd}
        if test $? -ne 0; then
          echo ${errmsg}
          echo "  python failed!!"
          echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn}                  \
            -f ${simsfcmapdiagsf} -v ${vadd}
          exit
          rm ${simsfcmapdiagsf}
        fi
      done
    fi
  fi
}

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: "$(echo ${timeS} | tr '!' ' ')
      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 Hovmoeller_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
    echo "    Sounding comparison on : "${sndstid}" by "${expl}
    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
      echo "    Time-Series comparison on : "${sti}" for "${sfcv}
      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
      echo "    Map comparison on : "$(echo ${timeS} | tr '!' ' ')" by "${expl}" for "${sfcv}
      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
}

function multi_WindRose_plot(){
# Function to plot a WindRose with data from multiple soundings
#   sti: id of the sounding station
#   stn: name of the sounding station to appear in plot ('!' for spaces)
#   ofigdir: folder for the output of the figure
#   files: files to use
#   labs: labels of the legend
#   lmark: marker of lines to use
#   timeS: time of the plot as it should appear in the figure ('!' for spaces)
#   timefS: time of the plot as it should appear in the figure's file
#   figheader: header of the figure's file name

    valuesfig=$1

    stid=`echo ${valuesfig} | tr '%' ' ' | awk '{print $1}'`
    stn=`echo ${valuesfig} | tr '%' ' ' | awk '{print $2}'`
    ofigdir=`echo ${valuesfig} | tr '%' ' ' | awk '{print $3}'`
    files=`echo ${valuesfig} | tr '%' ' ' | awk '{print $4}'`
    labs=`echo ${valuesfig} | tr '%' ' ' | awk '{print $5}'`
    lmark=`echo ${valuesfig} | tr '%' ' ' | awk '{print $6}'`
    timeS=`echo ${valuesfig} | tr '%' ' ' | awk '{print $7}'`
    timefS=`echo ${valuesfig} | tr '%' ' ' | awk '{print $8}'`
    figheader=`echo ${valuesfig} | tr '%' ' ' | awk '{print $9}'`

    Nfiles=`echo ${labs} | tr ',' ' ' | wc -w | awk '{print $1-1}'`
    frqm='1'
    iff=1
    while test ${iff} -le ${Nfiles}; do frqm=${frqm}'@10'; iff=`expr ${iff} + 1`; done

    values='linepoint;multicol;pres;-;'${lmark}';'${frqm}';spectral;auto:'${labs}':'
    values=${values}'WindRose!obs!,!sim!comparison!at!'${stid}'!on!'${timeS}'!UTC:'
    values=${values}${kfig}':auto:''True'

    ofign=${ofigdir}'/'${figheader}'_'${stid}'_'${timefS}'.'${kfig}
    if ${gscratch}; then rm ${ofign} >& /dev/null; fi
    if test ! -f ${ofign}; then
      echo "    WindRose comparison at: "${stid}" on "$(echo ${timeS} | tr '!' ' ')
      python ${pyHOME}/drawing.py -o draw_multiWindRose -f ${files} -S ${values}
      if test $? -ne 0; then
        echo ${errmsg}
        echo "  python failed!!"
        echo python ${pyHOME}/drawing.py -o draw_multiWindRose -f "'"${files}"'"     \
         -S "'"${values}"'"
        exit
      fi
      convert -trim multi_WindRose.${kfig} ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_multiWindRose -f "'"${files}"'"       \
        -S "'"${values}"'" >> ${ofilefigs}
      echo " " >> ${ofilefigs}
    fi
}

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

mkdir -p ${odir}
ofileins=${wdir}/allins_functions.inf
ofilefigs=${wdir}/allins_figures.inf
if ${fscratch}; then 
  echo "*** Instruction from all functions _______" > ${ofileins}
  echo " " >> ${ofileins}
fi
if ${gscratch}; then 
  echo "*** Instruction from all figures _______" > ${ofilefigs}
  echo " " >> ${ofilefigs}
fi

# Files in sounding folder (from UWyoming_netcdf.py)
if test ! ${sndHf} = 'None'; then
  sndfs=`ls -1 ${snddir}/${sndHf}* | tr '\n' ':'`
  echo "Soundings files _______"
  echo ${sndfs} | tr ':' '\n'
  sndfiles=`echo ${sndfs} | tr ':' ' '`
else
  sndfiles='None'
fi

# Files in sounding folder (from UWyoming_netcdf.py)
if test ! ${sfcHf} = 'None'; then
  sfcfs=`ls -1 ${sfcdir}/${sfcHf}** | tr '\n' ':'`
  echo "Surface files _______"
  echo ${sfcfs} | tr ':' '\n'
  sfcfiles=`echo ${sfcfs} | tr ':' ' '`
else
  sfcfiles='None'
fi

exps=`echo ${sims} | tr ':' ' '`
simH=`echo ${simHMT} | tr ',' ' ' | awk '{print $1}'`
simM=`echo ${simHMT} | tr ',' ' ' | awk '{print $2}'`
simT=`echo ${simHMT} | tr ',' ' ' | awk '{print $3}'`

###
## Preparing simulations
###

# sounding data
iost=1
obstype='sounding'
prn='snd'
if test ! ${sndHf} = 'None'; then 
  listfiles=${sndfiles}
else
  listfiles=' '
fi
wrfdiags=${wrfsnddiags}
simdiagd=${simsnddiagd}
for ostf in ${listfiles}; do
  stid=`python $pyHOME/nc_var.py -o grattr -f ${ostf} -S Station_number`
  ostlon=`ncdump -h ${ostf} | grep Station_longitude | awk '{print $3}'`
  ostlat=`ncdump -h ${ostf} | grep Station_latitude | awk '{print $3}'`
  ostheight=`ncdump -h ${ostf} | grep Station_elevation | awk '{print $3}'`

  if test ${iost} -eq 1; then 
    oststids=${stid}'@'${ostf}
    echo ${prn}","${stid}","${ostlon}","${ostlat}","${ostheight} > ${wdir}/stations.inf
  else 
    oststids=${oststids}':'${stid}'@'${ostf}
    echo ${prn}","${stid}","${ostlon}","${ostlat}","${ostheight} >> ${wdir}/stations.inf
  fi

  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}'`
    expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

    simfiles=`ls -1 ${simdir}/${expf}/${simH}*${simM}*${simT}`

    ofile=${odir}/simout_${prn}vars_${stid}_${expl}.nc
    if ${fscratch}; then rm ${ofile} >& /dev/null; fi
    if test ! -f ${ofile}; then
      for simf in ${simfiles}; do
        simfn=`basename ${simf}`
        fdate=`echo ${simfn} | tr '.' ' ' | awk '{print $1}' | tr '_' ' ' |          \
          awk '{print $3"_"$4}'`

        simsndptf=${odir}'/simout_vars_'${prn}'pt_'${stid}'_'${expl}'_'${fdate}'.nc'
        simsnddiagsf=${odir}'/simout_'${prn}'diags_'${stid}'_'${expl}'_'${fdate}'.nc'

        fvals=${ostlon}'#'${ostlat}'#'${simf}'#XLONG#XLAT#'${wrfdiags}'#'
        fvals=${fvals}${simdiagd}
        get_WRFsounding_data ${fvals}

      # end of files
      done

      # Concatenating files
      Nconfc=`ls -1 ${odir}/simout_${prn}diags_*${stid}*${expl}*.nc | wc -l |        \
        awk '{print $1}'`
      if test ${Nconfc} -gt 1; then 
        values=${odir}',Time,time'
        HMT='simout_'${prn}'diags_,'${stid}'_'${expn},'nc'
        python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values}      \
         -f ${HMT} -v all
        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 all
          exit
        fi
        echo " * Concatenating "${obstype}" data at point: "${ostlon}", "${ostlat} >> ${ofileins}
        echo python  ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S \
          "'"${values}"'" -f ${HMT} -v all >> ${ofileins}
        echo " " >> ${ofileins}
        mv netcdf_fold_concatenated_HMT.nc ${ofile}
      else
        ffilen=`ls -1 ${odir}/simout_${prn}diags_*${stid}*${expl}*.nc`
        cp ${ffilen} ${ofile}
      fi
    fi
    iexp=`expr ${iexp} + 1`
  # end of experiments
  done

  iost=`expr ${iost} + 1`
done
if test ! ${sndHf} = 'None'; then sndstids=${oststids}; fi

# surface data at station
iost=1
obstype='surface'
prn='sfc'
if test ! ${sfcHf} = 'None'; then 
  listfiles=${sndfiles}
else
  listfiles=' '
fi
wrfdiags=${wrfsfcdiags}
simdiagd=${simsfcdiagd}
for ostf in ${listfiles}; do
  stid=`python $pyHOME/nc_var.py -o varout -f ${ostf} -S 'time:0' -v id | grep NC \
    | awk '{printf("%d",$2)}'`
  ostlon=`python $pyHOME/nc_var.py -o varout -f ${ostf} -S 'time:0' -v lon | grep NC \
    | awk '{print $2}'`
  ostlat=`python $pyHOME/nc_var.py -o varout -f ${ostf} -S 'time:0' -v lat | grep NC \
    | awk '{print $2}'`
  ostheight=`python $pyHOME/nc_var.py -o varout -f ${ostf} -S 'time:0' -v height |   \
    grep NC | awk '{print $2}'`

  if test ! ${wdir}/stations.inf; then 
    oststids=${stid}'@'${ostf}
    echo ${prn}","${stid}","${ostlon}","${ostlat}","${ostheight} > ${wdir}/stations.inf
  else 
    oststids=${oststids}':'${stid}'@'${ostf}
    echo ${prn}","${stid}","${ostlon}","${ostlat}","${ostheight} >> ${wdir}/stations.inf
  fi

  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}'`
    expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

    simfiles=`ls -1 ${simdir}/${expf}/${simH}*${simM}*${simT}`

    for simf in ${simfiles}; do
      simfn=`basename ${simf}`
      fdate=`echo ${simfn} | tr '.' ' ' | awk '{print $1}' | tr '_' ' ' |            \
        awk '{print $3"_"$4}'`

      simsfcptf=${odir}'/simout_vars_'${prn}'pt_'${stid}'_'${expl}'_'${fdate}'.nc'
      simsfcdiagsf=${odir}'/simout_'${prn}'diags_'${stid}'_'${expl}'_'${fdate}'.nc'

      fvals=${ostlon}'#'${ostlat}'#'${simf}'#XLONG#XLAT#'${wrfdiags}'#'
      fvals=${fvals}${wrfsfcnondiags}'#'${simdiagd}
      get_WRFsurface_data ${fvals}

    # end of files
    done

    # Concatenating files
    ofile=${odir}/simout_${prn}vars_${stid}_${expl}.nc
    if ${fscratch}; then rm ${ofile} >& /dev/null; fi
    if test ! -f ${ofile}; then
      Nconfc=`ls -1 ${odir}/simout_${prn}diags_*${stid}*${expl}*.nc | wc -l | awk '{print $1}'`
      if test ${Nconfc} -gt 1; then 
        values=${odir}',Time,time'
        HMT='simout_'${prn}'diags_,'${stid}'_'${expn},'nc'
        python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values}      \
         -f ${HMT} -v all
        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 all
          exit
        fi
        echo " * Concatenating "${obstype}" data at point: "${ostlon}", "${ostlat} >> ${ofileins}
        echo python  ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S \
          "'"${values}"'" -f ${HMT} -v all >> ${ofileins}
        echo " " >> ${ofileins}
        mv netcdf_fold_concatenated_HMT.nc ${ofile}
      else
        ffilen=`ls -1 ${odir}/simout_${prn}diags_*${stid}*${expl}*.nc`
        cp ${ffilen} ${ofile}
      fi
    fi
    iexp=`expr ${iexp} + 1`
  # end of experiments
  done

  iost=`expr ${iost} + 1`
done
if test ! ${sfcHf} = 'None'; then sfcstids=${oststids}; fi

# surface data for maps
iost=1
obstype='surface_map'
prn='sfc'
listfiles=${sfcfiles}
wrfdiags=${wrfsfcdiags}
simdiagd=${simsfcmapdiagd}
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}'`
  expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

  simfiles=`ls -1 ${simdir}/${expf}/${simH}*${simM}*${simT}`

  for simf in ${simfiles}; do
    simfn=`basename ${simf}`
    fdate=`echo ${simfn} | tr '.' ' ' | awk '{print $1}' | tr '_' ' ' |            \
      awk '{print $3"_"$4}'`

    simsfcmapf=${odir}'/simout_vars_'${prn}'map_'${expl}'_'${fdate}'.nc'
    simsfcmapdiagsf=${odir}'/simout_'${prn}'mapdiags_'${expl}'_'${fdate}'.nc'

    fvals=${simf}'#'${wrfdiags}'#'${wrfsfcnondiags}'#'${simdiagd}
    get_WRFsfc2D_data ${fvals}

  # end of files
  done

  # Concatenating files
  ofile=${odir}/simout_${prn}mapvars_${expl}.nc
  if ${fscratch}; then rm ${ofile} >& /dev/null; fi
  if test ! -f ${ofile}; then
    Nconfc=`ls -1 ${odir}/simout_${prn}mapdiags_*${expl}*.nc | wc -l | awk '{print $1}'`
    if test ${Nconfc} -gt 1; then 
      values=${odir}',Time,time'
      HMT='simout_'${prn}'mapdiags_,'${expn},'nc'
      python ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S ${values}      \
       -f ${HMT} -v all
      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 all
        exit
      fi
      echo " * Concatenating "${obstype}" data: " >> ${ofileins}
      echo python  ${pyHOME}/nc_var.py -o netcdf_fold_concatenation_HMT -S           \
        "'"${values}"'" -f ${HMT} -v all >> ${ofileins}
      echo " " >> ${ofileins}
      mv netcdf_fold_concatenated_HMT.nc ${ofile}
    else
      ffilen=`ls -1 ${odir}/simout_${prn}mapdiags_*${expl}*.nc`
      cp ${ffilen} ${ofile}
    fi
  fi
  iexp=`expr ${iexp} + 1`
# end of experiments
done

###
##
# Figures
##
###

mkdir -p ${ofigdir}

##
# Stations maps
##
ofig=${ofigdir}/stations.png
if ${gscratch}; then rm ${ofig} >& /dev/null; fi
if test ! -f ${ofig}; then
  fivals=${mapstchar}':'${mapv}':lonlatxtrms:'${figexperimentname}'!stations:'
  fivals=${fivals}${kfig}':true'
  python ${pyHOME}/drawing.py -o draw_stations_map -f ${wdir}/stations.inf           \
    -S ${fivals}
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/drawing.py -o draw_stations_map -f ${wdir}/stations.inf    \
      -S "'"${fivals}"'"
    exit
  fi
  echo " * Map of the stations: " >> ${ofilefigs}
  echo python ${pyHOME}/drawing.py -o draw_stations_map -f ${wdir}/stations.inf      \
    -S "'"${fivals}"'" >> ${ofilefigs}
  echo " " >> ${ofilefigs}
  convert -trim stations_map.${kfig} ${ofig}
fi

###
## Soundings figures
###

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 'uavaFROMobswswd|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 'uavaFROMobswswd|ws@wd'
      exit
    fi
    echo " * Getting ua,va winds from station "${stid}": " >> ${ofileins}
    echo python ${pyHOME}/diagnostics.py -f ${sndf} -d "'"time@time,pres@pres"'"     \
      -v "'"uavaFROMobswswd|ws@wd"'" >> ${ofileins}
    echo " " >> ${ofileins}    
    mv diagnostics.nc ${ouavasndf}
  fi
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}"`

    sndstn=`python $pyHOME/nc_var.py -o grattr -f ${sndorigfn} -S Station_name |     \
      tr ' ' '!'`

    fivals=${it}'#'${stid}'#'${sndstn}'#'${timeS}'#'${ofigdir}'#'${timefS}'#'
    fivals=${fivals}${ouavasndf}'#WindRose_obs'
    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}`

    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 ',' ' '`
  sndstn=`python $pyHOME/nc_var.py -o grattr -f ${sndorigfn} -S Station_name |       \
    tr ' ' '!'`

  # 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}'`
      expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

      #simjoinselvars=${odir}'/simout_vars_'${expl}'.nc'
      #simdiagsf=${odir}'/simout_snddiags_'${expl}'.nc'
      #simfiles=`ls -1 ${simdir}/${expf}/${simH}*${simM}*${simT}`
      simjoinselvars=`ls -1 ${simdir}/${expf}/${simH}*${simM}*${simT} | head -n 1`
      simdiagsf=${odir}'/simout_sndvars_'${stid}'_'${expl}'.nc'
      simsndptf=${odir}'/simout_sndvars_'${stid}'_'${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

      # 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}'`
      timestepdiff=`echo ${timestepv} | tr '@' ' ' | awk '{print $2}'`
      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'

        # multi Wind-Rose
        # Getting pressure level from observations to plot
        if test ${it} -eq 0; then
          fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None        \
            -f ${simsndptf} -v pres | grep MAT | grep pres`
          nsimpres=`echo ${fstats} | awk '{printf ("%d",$3)}'`
        fi

        obsuavaf=${odir}/sounding_uava_${stid}.nc
        obspres=`python $pyHOME/nc_var.py -o varout -f ${obsuavaf} -S time:${it}     \
          -v pres | awk '{print $2}' | tr '\n' ':'`
        Lobspres1=`expr length ${obspres} | awk '{print $1-1}'`
        vals=${obspres:0:${Lobspres1}}',threshold,:,'${nsimpres}
        kobs=`python ${pyHOME}/generic.py -o inf_operSlist -S ${vals}`
        if test $? -ne 0; then
            echo ${errormsg}
            echo "  python failed !!"
            echo python ${pyHOME}/generic.py -o inf_operSlist -S ${vals}
            exit
        fi
        vals=${obspres:0:${Lobspres1}}',getvalue,:,'${kobs}
        onpres=`python ${pyHOME}/generic.py -o inf_operSlist -S ${vals}`
        echo "      obs level of sim npres "${nsimpres}": "${kobs}" "${onpres}

        mWRcomplabs='obs,'${expn}
        mWRcompfiles=${odir}'/sounding_uava_'${stid}'.nc@pres|1,'${kobs}',1;time|'
        mWRcompfiles=${mWRcompfiles}${it}'@ua,va#'${simsndptf}'@bottom_top|-1;Time|'
        mWRcompfiles=${mWRcompfiles}${timestep}'@ua,va'
        mWRmarks='.@'${expm}
      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'
        # multi Wind-Rose
        mWRcomplabs=${mWRcomplabs}','${expn}
        mWRcompfiles=${mWRcompfiles}'#'${simsndptf}'@bottom_top|-1;Time|'${timestep}
        mWRcompfiles=${mWRcompfiles}'@ua,va'
        mWRmarks=${mWRmarks}'@'${expm}
      fi

      if test ${it} -eq 0; then
        # 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'

          Hovmoeller_SkewT_map_plot ${fivals}

        # end sounding vars
        done
      fi

      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}'@'
    fivals=${fivals}'SkewT-logP_obs-sim_step'
    multi_SkewT_logP_plot ${fivals}

    # Wind Rose comparison
    fivals=${stid}'%'${sndstn}'%'${ofigdir}'%'${mWRcompfiles}'%'${mWRcomplabs}'%'
    fivals=${fivals}${mWRmarks}'%'${dateS}'%'${datefS}'%WindRose_obs-sim_step'

    multi_WindRose_plot ${fivals}

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

# End of soundings stations
done

# Joining sounding 1D variables
obssnd1Dfile=${odir}'/all_sounding_1D.nc'
if test ! -f ${obssnd1Dfile}; then
  python ${pyHOME}/nc_var.py -o join_sounding_obsfiles -S ${snddir}':UWyoming_snd'   \
    -v all
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/nc_var.py -o join_sounding_obsfiles                        \
      -S ${snddir}':UWyoming_snd' -v all
  fi
  echo " * Joining all 1D variables from all sounding station files : " >> ${ofileins}
  echo python ${pyHOME}/nc_var.py -o join_sounding_obsfiles                          \
      -S ${snddir}':UWyoming_snd' -v all >> ${ofileins}    
  echo " " >> ${ofileins}    
  mv joined_soundingstations.nc ${obssnd1Dfile}
fi

exit

###
## SFC figures
###

# Joining all observations
obssfcfile=${odir}'/all_single-stations.nc'
if test ! -f ${obssfcfile}; then
  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles -S                       \
    ${sfcdir}':OBSnetcdf' -v all
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo  python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles \
      -S ${sfcdir}':OBSnetcdf' -v all
  fi
  echo " * Joining all surface station files : " >> ${ofileins}
  echo python ${pyHOME}/nc_var.py -o join_singlestation_obsfiles -S                  \
    ${sfcdir}':OBSnetcdf' -v all >> ${ofileins}    
  echo " " >> ${ofileins}    
  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}'`
      for stidf in ${sndids}; do
        stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
        simsfcdiagf=${odir}'/simout_sfcvars_'${stdi}'_'${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
        iexp=`expr ${iexp} + 1`
      done
    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}'`
      for stidf in ${sndids}; do
        stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
        simsfcdiagf=${odir}'/simout_sfcvars_'${stdi}'_'${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
    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}'`      
      for stidf in ${sndids}; do
        stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
        simsfcdiagf=${odir}'/simout_sfcvars_'${stdi}'_'${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
        iexp=`expr ${iexp} + 1`
      done
    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}'`
      for stidf in ${sndids}; do
        stid=`echo ${stidf} | tr '@' ' ' | awk '{print $1}'`
        simsfcdiagf=${odir}'/simout_sfcvars_'${stdi}'_'${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
    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}'`
    expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

    simsfcdiagf=${odir}'/simout_sfcvars_'${stdi}'_'${expl}'.nc'

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

    sfccompfiles=${sfccompfiles}','${simsfcdiagf}'%time|-1'
    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}
    exit

  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}'`
  expm=`echo ${exp} | tr ',' ' ' | awk '{print $5}'`

  simsfcdiagf=${odir}'/simout_sfcmapdiags_'${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