#!/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
#

# files from scratch
fscratch=false

# figures from scratch
gscratch=true

# Working directory
wdir='/home/lluis/estudios/ChemGBsAs/tests/199501'

## U. Wyoming file (respect wdir)
snddir='obs/snd'
sndorigfn='UWyoming_snd_87576.nc'
snddimt=6
sndstn='Ezeiza'

# Kind of graphical file output
kfig='png'

# Color for SkewT_evol
col=255
col2=200
fraccol=16
fraccol2=16

# Simulation dir
simdir='/home/lluis/PY'

# Variables from the simulations
simvars='XLONG,XLAT,T,P,PB,PH,PHB,QVAPOR,U,V,MU,MUB,SINALPHA,COSALPHA,T2,U10,V10,'
simvars=${simvars}'PSFC,Q2'

# Diagnotics for the soundigns
wrfsnddiags='WRFt:WRFtda:WRFp:WRFua:WRFva:ws'

# Figures output directory
ofigdir=${wdir}/figs

# Sounding comparison vars
sndvars='ta:tda:ws'

# Observations surface files
sfcobsdir='obs/sfc_CAM/'
obsHf='OBSnetcdf_'
# Surface variables [varn1]@[minv1]@[maxv1]:[varn2]@[minv2]@[maxv2]:...
#sfcvars='tas:wds:wss:psl:clt:tdas'
sfcvars='tas@FROMobs@FROMobs:wds@0.@6.28:wss@0.@@FROMobs:tdas@FROMobs@FROMobs'
#wrfsfcdiags='wds:wss:WRFpsl_ecmwf:clt:WRFtdas'
wrfsfcdiags='T2:wds:wss:WRFtdas'

# map values
mapv='cyl,f'

# map coverage for cont-disc
mapcover='sponge,0.1,0.1'

# Format of time-axis for time-series
fmtTts='exct,12,h|%d$^{%H}$'

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

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

#######    #######
## MAIN
    #######
ofilefigs=${wdir}/allins_figures.inf
echo "*** INstruction from all figures" > ${ofilefigs}

mkdir -p ${ofigdir}

# Soundings plots
cd ${wdir}/${snddir}

# 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}

# Sounding

# Getting wind components
ofile='sounding_uava.nc'
if ${fscratch}; then rm ${ofile}; fi
if test ! -f ${ofile}; then
  python ${pyHOME}/diagnostics.py -f ${sndorigfn} -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 ${sndorigfn} -d 'time@time,pres@pres' -v 'uavaFROMwswd|ws@wd'
    exit
  fi 
  mv diagnostics.nc ${ofile}
fi

it=0
tatda_evol_values=''
tatda_evol_files=''
sndtimes=''
snddates=''
sndFdates=''
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 +%Y"/"%m"/"%d"!"%H -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}"`
  # Wind rose
  values='pres|-1;time|'${it}':linepoint;multicol;pres;auto;auto;rainbow;auto:'
  values=${values}${sndstn}'!sounding!WindRose!on!'${timeS}'!local!time:png:cardinals:False:WindRose:True'
  ofign=${ofigdir}'/WindRose_'${timefS}'.png'
  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 ${ofile} >& /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 ${ofile}
      exit
    fi
    mv WindRose.png ${ofign}
    echo "* "${ofign} >> ${ofilefigs}
    echo python ${pyHOME}/drawing.py -o draw_WindRose -S ${values} -v ua,va -f ${ofile} >> ${ofilefigs}
  fi

  # Sounding
  echo ${sndorigfn}
  values='time|'${it}',pres|-1:auto:auto:'${sndstn}'!sounding!on!'${timeS}'!local!time:png:True'
  ofign=${ofigdir}'/SkewT-logP_'${timefS}'.png'
  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 ${sndorigfn} >& /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 ${sndorigfn}
      exit
    fi
    mv SkewT.png ${ofign}
    echo "* "${ofign} >> ${ofilefigs}
    echo python ${pyHOME}/drawing.py -o draw_SkewT -S ${values} -v ta,tda,pres -f ${sndorigfn} >> ${ofilefigs}
  fi

  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} - ${fraccol}`
  col2=`expr ${col2} - ${fraccol2}`

  #exit

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

# multi_sounding
tatda_evol_values='auto:auto:multilines!'${tatda_evol_labs}'!'${tatda_evol_cols}'!'
tatda_evol_values=${tatda_evol_values}'-!,!2:0,auto:!sounding!evolution!on!'
tatda_evol_values=${tatda_evol_values}${timeS3}'!local!time:png:True'
ofign=${ofigdir}'/SkewT-logP_evol.png'
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!!"
    python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${tatda_evol_values} -f ${tatda_evol_files}
    exit
  fi
  mv 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}
fi

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

# Joining files
ofile0='simout_vars.nc'
if test ! -f ${ofile0}; then
  values=${simdir}',Time,WRFtime'
  HMT='wrfout,'${timeS4},'00'
  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 ${wdir}/${ofile0}
fi

gridsnd=`python ${pyHOME}/nc_var.py -o get_point -f ${wdir}/${ofile0}                \
  -S 'XLONG:XLAT:Time|0' -v ${sndlon}','${sndlat}`
xsnd=`echo ${gridsnd} | tr ',' ' ' | awk '{print $1}'`
ysnd=`echo ${gridsnd} | tr ',' ' ' | awk '{print $2}'`
echo "Equivalent simulated sounding grid point: "${gridsnd}

# Computing sounding diagnostics
diagns=`echo ${wrfsnddiags} | tr ':' ' '`
diagd='time@time,bottom_top@P,south_north@XLAT,west_east@XLONG'
diagvals=''
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
    diagvals=${diagv}
  else
    diagvals=${diagvals}','${diagv}
  fi
  idiag=`expr ${idiag} + 1`
done

ofile='simout_snddiags.nc'
if test ! -f ${ofile}; then
  python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
    exit
  fi
  mv diagnostics.nc simout_snddiags.nc
  #To deg
  python ${pyHOME}/nc_var.py -o valmod -S subc,273.15 -f ${ofile} -v ta
  #To Pa
  python ${pyHOME}/nc_var.py -o valmod -S divc,100. -f ${ofile} -v pres
fi

# Getting values at the sounding point
ofile2='simout_vars_sndpt.nc'
if test ! -f ${ofile2}; 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 ${wdir}/${ofile}         \
    -S ${values} -v all
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/nc_var.py -o DataSetSection_multidims -f ${wdir}/${ofile}  \
    -S ${values} -v all
    exit
  fi
  ofilen=`ls -rt1 | tail -n 1`
  mv ${ofilen} ${ofile2}
fi

# Time information from simulation
simdimt=`python ${pyHOME}/nc_var.py -o itime -S CFtime, -v time -f simout_snddiags.nc\
  | grep dimt | awk '{print $2}'`
simtunits=`python ${pyHOME}/nc_var.py -o ivattrs -v time -f simout_snddiags.nc |     \
  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}

# Plotting
complabs='$ta^{obs}$,$ta^{sim}$,$tda^{obs}$,$tda^{sim}$'
compcols='#FF6464,#FF6464,#6464FF,#6464FF'
it=0
for its in ${sndts}; do
  timestep=`python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${ofile2}              \
    -S ${its}';'${tunits} -v time`
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${ofile2}                 \
      -S ${its}';'${tunits} -v time
  fi

  dateS=`echo ${snddates} | tr ',' '\n' | head -n ${it} | tail -n 1`
  datefS=`echo ${sndFdates} | tr ',' '\n' | head -n ${it} | tail -n 1`
  echo ${its}": it="${timestep}" "$(echo ${dateS} | tr '!' ' ')" "${datefS}

  values='auto:auto:multilines!'${complabs}'!'${compcols}'!'
  values=${values}'-,-.,-,-.!,!2:0,auto:!sounding!evolution!on!'
  values=${values}${dateS}'!local!time:png:True'

  compfiles=${wdir}/${snddir}/${sndorigfn}':pres|-1,time|'${it}':ta,pres;'
  compfiles=${compfiles}${ofile2}':bottom_top|-1,time|'${timestep}':ta,pres;'
  compfiles=${compfiles}${wdir}/${snddir}/${sndorigfn}':pres|-1,time|'${it}':tda,pres;'
  compfiles=${compfiles}${ofile2}':bottom_top|-1,time|'${timestep}':tda,pres'

  ofign=${ofigdir}'/SkewT-logP_obs-sim_'${datefS}'.png'
  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
    mv multi_SkewT.png ${ofign}
    echo "* "${ofign} >> ${ofilefigs}
    echo python ${pyHOME}/drawing.py -o draw_multi_SkewT -S ${values} -f ${compfiles} >> ${ofilefigs}

  fi
  it=`expr ${it} + 1`
done

# Plotting snd comparison

# Soundings
sndvs=`echo ${sndvars} | tr ':' ' '`
for sndv in ${sndvs}; do
  ofign=${ofigdir}'/SkewT-logP_obs-sim_evol_'${sndv}'.png'

  cfiles='simout_vars_sndpt.nc;'${sndv}';time;pres;time|-1,pres|-1@'
  cfiles=${cfiles}'obs/snd/UWyoming_snd_87576.nc;'${sndv}';time;pres;time|-1,pres|-1'
  cvalues=${sndv}';y;auto|exct,12,h|%d$^{%H}$|time!($[DD]^{[HH]}$);Vfix,auto,50.,auto'
  cvalues=${cvalues}';auto;Srange,Srange;auto;obs!&!sim!Ezeiza!airport!sounding;png;'
  cvalues=${cvalues}'flip@y;None;yes'
  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
    mv 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}
  fi
done

# Plotting map comparison
cd $wdir

# Joining all observations
obssfcfile='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 ${wdir}/${sfcobsdir}/${obssfcfile}
fi

# Computing surface diagnostics
diagns=`echo ${wrfsfcdiags} | tr ':' ' '`
diagd='time@time,south_north@XLAT,west_east@XLONG'
diagvals=''
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
      diagvals=${diagv}
      idiag=`expr ${idiag} + 1`
    else
      diagvals=${diagvals}','${diagv}
    fi
  fi
done

osfcdiagfile='simout_sfcdiags.nc'
if test ! -f ${osfcdiagfile}; then
  python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
  if test $? -ne 0; then
    echo ${errmsg}
    echo "  python failed!!"
    echo python ${pyHOME}/diagnostics.py -f ${ofile0} -d ${diagd} -v ${diagvals}
    exit
  fi
  mv diagnostics.nc ${osfcdiagfile}
  #To rad
  python ${pyHOME}/nc_var.py -o valmod -S mulc,0.0174532925199 -f ${osfcdiagfile} -v wds
  python ${pyHOME}/nc_var.py -o varaddattr -S 'units|rad' -f ${osfcdiagfile} -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 ${ofile0},${vadd} -f ${osfcdiagfile}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o fvaradd -S ${osfcdiagfile},${vadd} -f ${ofile0}
      exit
      rm ${osfcdiagfile}
    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 ${osfcdiagfile}
    fi
    python $pyHOME/nc_var.py -o chvarname -S ${CFvarn} -f ${osfcdiagfile} -v ${vadd}
    if test $? -ne 0; then
      echo ${errmsg}
      echo "  python failed!!"
      echo python ${pyHOME}/nc_var.py -o chvarname -S ${CFvarn} -f ${osfcdiagfile} -v ${vadd}
      exit
      rm ${osfcdiagfile}
    fi
  done
fi

sfcs=`echo ${sfcvars} | tr ':' ' '`

# getting min,max values
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 ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc=`echo ${fstats} | awk '{print $3}'`
  elif test ${nsfc} = 'FROMsims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc=`echo ${fstats} | awk '{print $3}'`
  elif test ${nsfc} = 'FROMobssims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc1=`echo ${fstats} | awk '{print $3}'`
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc2=`echo ${fstats} | awk '{print $3}'`
    echo ${nsfc1}" "${nsfc2} | awk '{if ($1 < $2) {print $1;} else {print $2;}}'
  fi
  if test ${xsfc} = 'FROMobs'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    xsfc=`echo ${fstats} | awk '{print $4}'`
  elif test ${xsfc} = 'FROMsims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    xsfc=`echo ${fstats} | awk '{print $4}'`
  elif test ${nsfc} = 'FROMobssims'; then
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${wdir}/${sfcobsdir}/${obssfcfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc1=`echo ${fstats} | awk '{print $3}'`
    fstats=`python $pyHOME/nc_var.py -o field_stats -S full,1.e+20,None              \
      -f ${osfcdiagfile} -v ${sfcv} | grep MAT | grep ${sfcv}`
    nsfc2=`echo ${fstats} | awk '{print $3}'`
    echo ${nsfc1}" "${nsfc2} | awk '{if ($1 > $2) {print $1;} else {print $2;}}'
  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 \
  ${wdir}/${sfcobsdir}/${obssfcfile} -S 'Nstations:-1|time:0' | \
  awk '{$1=""; printf("%d:",substr($0,2))}'`
stsnames=`python ${pyHOME}/nc_var.py -o varout -v stsname -f \
  ${wdir}/${sfcobsdir}/${obssfcfile} -S Nstations:-1 | \
  awk '{$1=""; print substr($0,2)}' | tr '\n' ':' | tr ' ' '!'`
stslons=`python ${pyHOME}/nc_var.py -o varout -v stslon -f \
  ${wdir}/${sfcobsdir}/${obssfcfile} -S Nstations:-1 | \
  awk '{$1=""; print substr($0,2)}' | tr '\n' ':'`
stslats=`python ${pyHOME}/nc_var.py -o varout -v stslat -f \
  ${wdir}/${sfcobsdir}/${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`
  gridsfc=`python ${pyHOME}/nc_var.py -o get_point -f ${osfcdiagfile}              \
    -S 'XLONG:XLAT:Time|0' -v ${stl}','${stL}`
  xgrid=`echo ${gridsfc} | tr ',' ' ' | awk '{print $1}'`
  ygrid=`echo ${gridsfc} | tr ',' ' ' | awk '{print $2}'`
  sfcstats=${sti}'@'${stn}'@'${stl}'@'${stL}'@'${gridsfc}
  ist_1=`expr ${ist} - 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}'`

    ofign=${ofigdir}'/obs-sim_evol_'${sfcv}_${sti}'.'${kfig}
    files=${wdir}/${sfcobsdir}/${obssfcfile}'%time|-1;Nstations|'${ist_1}','
    files=${files}${osfcdiagfile}'%time|-1;south_north|'${ygrid}';west_east|'${xgrid}
    values='ststimes,time;y;time!([DD]$^{[HH]}$);auto;obs,sim;'${sfcv}';'${sfcv}
    values=${values}'|evolution|at|'$(echo ${stn} | tr '!' '|')';'${nsfc}','${xsfc}';'
    values=${values}'time|'${fmtTts}';0|12;'
    values=${values}${kfig}';-;k,b;.;2.;2.;all;-1;True'
    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!!"
        python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv}
        exit
      fi
      mv lines_time.${kfig} ${ofign}
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_lines_time -f ${files} -S ${values} -v ${sfcv} >> ${ofilefigs}
    fi
    #exit
  done
  ist=`expr $ist + 1`
  #exit
done
exit

# Cont-disc comparison
it=0
while test ${it} -le ${simdimt}; do
  simtstep=`python ${pyHOME}/nc_var.py -o varout -f ${wdir}/${osfcdiagfile} -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 ${wdir}/${osfcdiagfile} -S          \
      'time:'${it} -v time
  fi

  timeS=`date +%Y"/"%m"/"%d"!"%H -d"${simDref} ${simtstep} ${simTu}"`
  timeS2=`date +%d"$^{"%H"}$" -d"${simDref} ${simtstep} ${simTu}"`
  timeS3=`date +%Y"-"%b -d"${simDref} ${simtstep} ${simTu}"`
  timeS4=`date +%Y -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
  obsf=${wdir}/${sfcobsdir}/${obsfile}
  oit=`python ${pyHOME}/nc_var.py -o get_time -f ${wdir}/${sfcobsdir}/${obssfcfile}  \
    -S ${simtstep}';'${tunits} -v ststimes | grep -v differ`

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

    ofign=${ofigdir}'/map_obs-sim_'${timefS}'_'${sfcv}'.png'

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

    cfiles=${wdir}/${osfcdiagfile}';'${sfcv}';XLONG;XLAT;Time|'${it}',time|'${it}','
    cfiles=${cfiles}'west_east|-1,south_north|-1@'${wdir}/${sfcobsdir}/${obssfcfile}
    cfiles=${cfiles}';'${sfcv}';stslon;stslat;time|'${oit}',lon|-1,lat|-1'
    cvalues=${sfcv}':west_east,south_north:auto:'${cbar}',auto,auto:'${nsfc}','
    cvalues=${cvalues}${xsfc}':''auto:obs!&!sim!'${sfcv}'!on!'${timeS}'!UTC:png:None:'
    cvalues=${cvalues}${mapv}':'${mapcover}':yes'

    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
      mv 2Dshad_obs-sim_comparison.png ${ofign} 
      echo "* "${ofign} >> ${ofilefigs}
      echo python ${pyHOME}/drawing.py -o draw_2D_shad_contdisc -f ${cfiles} -S ${cvalues} >> ${ofilefigs}

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