#!/bin/bash

########################################################################
# Compute average btwn t1,t2,zmin,zmax for ncdf file
# Auteur: C Rio, F Hourdin
# Modified: F Couvreux
# Modified: 23/01/2019 N Villefranque from r110
########################################################################

set +vx 
metric_full_name=$1


dirWAVE=$2

# Erreur sur la température
dtmin=0.1 # L95
dtmin=0.3

# Erreurs sur l'humidité spécifique
dqmin=0.0005 # L95
dqmin=0.0007

# Erreur relative sur les couvertures nuageuses (moyennes ou max)
err_neb=0.3 # L95
err_neb=0.3

# Erreur sur les hauteurs de nuages (a relier a la resolution verticale du SCM)
fact_dz_min=0.1 # L95
fact_dz_min=0.2

#########################################################################
# Decompisition of inputs
#########################################################################
IFS='_' read -a temp <<< $metric_full_name
# Decomposition soit en
# cas_souscas_metrique_t1_t2 => comparaison SCM / LES ou LES / LES
# RAD_cas_souscas_metrique_t2_t2 => comparaison ecRad / MC sur le champ de cas/souscas, un souscas par instant
# Attention de ne pas appeler un cas "RAD"
tmp=${temp[0]}
if [ ${tmp} == RAD ] # 3 first characters
then
  REF=RAD
  case_basis=${temp[1]}
  name_subcase=${temp[2]}
  metric_name=${temp[3]}
  tmin=${temp[4]}
  tmax=${temp[5]} 
  INS=${REF}${name_subcase: -3}
else
  REF=LES
  case_basis=${temp[0]}
  name_subcase=${temp[1]}
  metric_name=${temp[2]}
  tmin=${temp[3]}
  tmax=${temp[4]} 
fi

# metric_output_name : LES_cas_souscas_... ou RAD_cas_souscas_... ou WAVE1_cas_souscas...
metric_output_name=${dirWAVE}_$metric_full_name
casename=$case_basis/$name_subcase
echo Cas $casename, metric $metric_name, time $tmin : $tmax, treating $dirWAVE
echo $metric_output_name # same info as above

case $case_basis in
	ARMCU|RICO|IHOP) zmax_domain=3500 ;;
	SANDU) zmax_domain=2800 ;;
	*) zmax_domain=2000
esac


#########################################################################
nWAVE=${dirWAVE:4}

# a changer en fonction du repertoire de l'utilisateur
DIRin=$dirWAVE/$casename

########################################################################
# list of files to be post processed
# LES or SCM
# trick to compute LES std when one LES only is available
# by shifting the time of the metrics by +1 or -1
# Default : no shift (shf=+0)
########################################################################
cd $DIRin
shfs="+0"
if [ -f LES0.nc ] ; then
   prefile=LES
   if [ ! -f LES1.nc ] ; then shfs="+0 -1 +1" ; fi
elif [ -f LESLESSCM_${nWAVE}-001.nc ] ; then
   prefile=LESLESSCM_${nWAVE}-
elif [ -f ${INS}.nc ] ; then 
   prefile=RAD
else
   prefile=SCM_${nWAVE}-
fi

ls ${prefile}*nc >list
sed -e "s/$prefile//" -e "s/.nc//" list > list_nruns
nruns=`wc -l list_nruns | awk '{print $1}'`
nr=0
cd -

########################################################################
# Loop on metrics
########################################################################
tmp=metrics_tmp.csv
echo $metric_output_name >| $tmp
for run in `more $DIRin/list_nruns`; do
  for shf in $shfs ; do
   run3=${run}
   file=$DIRin/$prefile${run3}
   (( nr = $nr + 1 ))
   #################################################################
   #defini les bornes temporelles et verticales
   #################################################################
   t1=`echo $tmin | awk ' { print $1 - 1 '$shf' } '`
   t2=`echo $tmax | awk ' { print $1 - 1 '$shf' } '`
   t1R=`echo $tmin | awk ' { print $1 '$shf' } '`
   t2R=`echo $tmax | awk ' { print $1 '$shf' } '`
   case ${metric_name:0:3} in

      #################################################################
      #Metrics computed in R with htune_netcdf2csvMetrics.R
      "neb"|"lwp"|"Ay-"|"net"|"rat"|"dns"|"tra"|"upt") 
      #################################################################
	  metric=`Rscript --vanilla htune_netcdf2csvMetrics.R $file.nc $metric_name $t1R $t2R $zmax_domain | awk ' { print $2 } '`
	  ;;

      #################################################################
      #Metrics computed with cdo
      "zav")
      #################################################################
	  z1=`echo $metric_name | awk -F- ' { print $2 } '`
	  z2=`echo $metric_name | awk -F- ' { print $3 } '`
	  var=`echo $metric_name | awk -F- ' { print $4} '`
	   
          rm -f XXXX_intermediate_extract*.nc
          timename=`ncdump -h $file.nc | grep -i time |head -1 | awk ' {print $1 } '`
          ncap2 -s "mask= (zf >= ${z1} && zf <= ${z2})" ${file}.nc XXXX_intermediate_extract.nc
          #echo "execute ncwa to perform time and z average"
          ncwa -m mask -a ${timename} -d ${timename},${t1},${t2} -v ${var} XXXX_intermediate_extract.nc XXXX_intermediate_extract2.nc
          #echo "execute ncwa to get the only non zero value"
          ncwa -y max XXXX_intermediate_extract2.nc XXXX_intermediate_extract3.nc
          # "put the value in type_metric_value"
          metric=`ncks -v ${var} -c -H XXXX_intermediate_extract3.nc | grep ${var} | grep '=' |  awk ' {print $3}'`
          ;;

      #################################################################
      *) echo Metrics $metric_name not available yet. ; echo 'Want to contribute ?' ; exit 1
      #################################################################

   esac
   echo $metric >> $tmp
   echo $metric
  done
done

####################################################################
# For LES : mean value and error
####################################################################

if [ $prefile = $REF ] ; then
  head=`head -1 $tmp`
  if [ $REF = LES ] ; then
    mean=`sed -n 2p $tmp`
    echo MEAN $mean
    std=`sed -n -e '2,$p' $tmp | awk ' BEGIN { std = 0 ; n=0 } { p = $1 - '$mean' ; std = std + p * p ; n = n + 1 } END { print std / ( n - 1 ) } '`
    case ${metric_name:0:3}  in
      # On impose que l'erreur sur la hauteur soit au moins 0.1 * Zmoy
      "neb") if [ "${metric_name}" = "nebmax" -o "${metric_name}" = "nebdz"  ] ; then fact=$err_neb ; else fact=$fact_dz_min ; fi ; std=`echo $std $mean $fact | awk ' { min=( $2 * $3 ) ^2 ; if ( $1 > min )  { print $1 } else { print min } } '` ;;
      "zav") case $var  in
        # On impose que l'erreur sur theta soit au moins 0.1 * THmoy et sur qv 0.0005*qvmoy
        "theta") std=`echo $std $mean | awk ' { min=( '$dtmin' ) ^2 ; if ( $1 > min )  { print $1 } else { print min } } '` ;;
        "qv") std=`echo $std $mean | awk ' { min=( '$dqmin' ) ^2 ; if ( $1 > min )  { print $1 } else { print min } } '` ;;
        *) ;;
      *) ;;
    esac
    esac
  else
    nbt=10#${INS:3} # need to force base 10 otherwise bash takes it as octal number
    line=$(( $nbt+1 )) # line is the timestep + 1 (the header)
    #mean=`sed -n ${line}p $tmp`
    mean=`sed -n 2p $tmp`
    file=${REF}/${casename}/${INS}
    std=`Rscript --vanilla htune_netcdf2csvMetrics.R $file.nc std_$metric_name $t1R $t2R $zmax_domain | awk ' { print $2 } '`
  fi
  echo TYPE,$head >| $metric_output_name.csv
  echo MEAN,$mean >> $metric_output_name.csv
  echo VAR,$std >> $metric_output_name.csv
else
  echo SIM > $DIRin/list_simus
  sed -e "s/.nc//" $DIRin/list >> $DIRin/list_simus
  paste -d, $DIRin/list_simus $tmp > $metric_output_name.csv
fi

\cp -f $metric_output_name.csv $DIRin

rm -f XXXX*.nc
rm -f *_metric_value
rm -f list*