#!/bin/sh
#
#****s* fm/test_1dvar_iono.sh
#
# NAME
#   test_1dvar_iono.sh - Test the ropp_1dvar program.
#
# SYNOPSIS
#   test_1dvar_iono.sh
#
# DESCRIPTION
#   The background is passed through ropp_fm -direct_ion.  The resulting
#   pseudo observations {bangle_L1, bangle_L2} are perturbed a bit.
#   The background {T, q, p*, Ne_max, H_peak, H_width} are also perturbed a bit.
#   A ropp-1dvar retrieval is then made on these two files.  The resulting
#   ionospheric parameters are compared to the original ones (for interest)
#   and to reference ones (for validation).
#
# USES
#   ropp_1dvar_bangle
#   (ncap2 to perturb the files - not needed by the general user)
#
# NOTES
#   1) Assumes PWD = ropp_1dvar/tests
#
# REFERENCES
#   See the ROPP 1DVAR User Guide (SAF/ROM/METO/UG/ROPP/007)
#
#****

if [ "$(echo $(basename $(dirname $PWD)) |cut -c1-10)"/"$(basename $PWD)" != "ropp_1dvar/tests" ] ; then
  echo "*** Not in ropp_1dvar*/tests subdirectory - test NOT PERFORMED"
  exit
fi

ulimit -S -s unlimited  # This test needs a stack size over 8192 Kibytes.

TEST=t_1dvar_iono_bangle

COMMENT="1DVAR L1 and L2"

EXTRA_CMD=" -direct_ion"

LOGFILE=${TEST}.log

echo " "
echo "Running $TEST ($COMMENT) ..."

echo "*** Results log of $TEST ($COMMENT) ***" > $LOGFILE


#-------------------------------------------------------------------------------
# 0. Loop over tests
#-------------------------------------------------------------------------------

for BGFILE in ../../ropp_fm/data/bgr20090401_000329_M02_2030337800_N0007_YYYY.nc ; do


# 1. Forward model to generate pseudo obs
# ---------------------------------------

# This is how to generate the dataset; for make tests it's stored in ../data

#  EXEC=../../ropp_fm/tools/ropp_fm_bg2ro_1d

#  OBFILE=$(basename $BGFILE |sed -es/'bgr'/'obs'/)

#  $EXEC  -direct_ion  -f  -d  $BGFILE  -o $OBFILE

#  cp $OBFILE ../data/$(basename $BGFILE |sed -es/'bgr'/'obs'/)

  OBFILE=../data/$(basename $BGFILE |sed -es/'bgr'/'obs'/)


# 2. Perturb the background {T, q, p*, Ne_max, H_peak, H_width}
# -------------------------------------------------------------

# This is how to generate the dataset; for make tests it's stored in ../data

#  BGFILE1=$(basename $BGFILE |sed -es/'.nc'/'1.nc'/)
#  ncap2 -h -s"temp=float(trunc(temp)); shum=float(trunc(shum)); press=float(1+10*trunc(press*0.1)); \
#          Ne_max=float(Ne_max+1.e10); H_peak=float(H_peak-20000); H_width=float(H_width+5000); \
#          bg_year=0*bg_year+2009; bg_month=0*bg_month+4; bg_day=0*bg_day+1; bg_hour=0*bg_hour+0; bg_minute=0*bg_minute+0" \
#          $BGFILE -O $BGFILE1a
#  ncap2 -h -s"where(shum < 0.5) shum=float(1.e-6)" \
#          $BGFILE1a -O $BGFILE1

#  cp $BGFILE1 ../data/$(basename $BGFILE |sed -es/'.nc'/'1_reference.nc'/)

  BGFILE1=../data/$(basename $BGFILE |sed -es/'.nc'/'1_reference.nc'/)


# 3. Perturb the obs {bangle_L1, bangle_L2}
# -----------------------------------------

# This is how to generate the dataset; for make tests it's stored in ../data

#  OBFILE1=$(basename $OBFILE |sed -es/'.nc'/'1.nc'/)

#  ncap2 -s"bangle_L1=1.e-5*trunc(1.e5*bangle_L1); \
#           bangle_L2=1.e-5*trunc(1.e5*bangle_L2)" $OBFILE -O $OBFILE1

#  cp $OBFILE1 ../data/$(basename $OBFILE |sed -es/'.nc'/'1_reference.nc'/)

  OBFILE1=../data/$(basename $OBFILE |sed -es/'.nc'/'1_reference.nc'/)


# 4. Carry out retrieval using this background. Can the minimiser find the original state?
# ----------------------------------------------------------------------------------------

  EXEC=../tools/ropp_1dvar_bangle

  BGCOV=../errors/ropp_bg_ecmwf_error_corr_L91.nc

  OBCOV=../errors/ropp_ob_bangle_error_corr_300L.nc

  CONFIG=../config/ecmwf_bangle_1dvar_iono.cf  # Uses bg_covar_method=FSFC, obs_covar_method=FSFC

  OFILE=$(basename $BGFILE |sed -es/'bgr'/'anl'/ |sed -es/'.nc'/'_iono.nc'/)

  if [ -f $EXEC ]; then

    echo "./$EXEC  ${EXTRA_CMD}  -d  -y $OBFILE1  --obs-corr $OBCOV  -b $BGFILE1  --bg-corr $BGCOV  -c $CONFIG  -o $OFILE" >> $LOGFILE
          ./$EXEC  ${EXTRA_CMD}  -d  -y $OBFILE1  --obs-corr $OBCOV  -b $BGFILE1  --bg-corr $BGCOV  -c $CONFIG  -o $OFILE >> $LOGFILE

  else

    echo "*** $EXEC not found - test FAILED"  >> $LOGFILE
    exit 1

  fi


# 5. Try a standard neutral retrieval on the same files
# -----------------------------------------------------

#  OFILE2=$(basename $BGFILE |sed -es/'bgr'/'anl'/ |sed -es/'.nc'/'_neut.nc'/)

#  echo " "
#  echo "$EXEC               -d  -y $OBFILE1  --obs-corr $OBCOV  -b $BGFILE1  --bg-corr $BGCOV  -c $CONFIG  -o $OFILE2"
#        $EXEC               -d  -y $OBFILE1  --obs-corr $OBCOV  -b $BGFILE1  --bg-corr $BGCOV  -c $CONFIG  -o $OFILE2


# 6. Summarise results
# --------------------

#  echo "Initially:"                         
#  ncks -Q -H -vNe_max,H_peak,H_width $BGFILE

#  echo "After perturbation:"                 
#  ncks -Q -H -vNe_max,H_peak,H_width $BGFILE1

#  echo "After ionospheric retrieval:"      
#  ncks -Q -H -vNe_max,H_peak,H_width $OFILE
##  ncks -Q -H -vJ,J_init,J_scaled $OFILE   

##  echo "After neutral retrieval:"           
##  ncks -Q -H -vNe_max,H_peak,H_width $OFILE2
##  ncks -Q -H -vJ,J_init,J_scaled $OFILE2    


# 7. Compare retrieved iono params with reference ones
# ----------------------------------------------------

  OFILE_REF=../data/$(basename $OFILE |sed -es/'.nc'/'_reference.nc'/)
  
# 7.1 Check output files produced
# -------------------------------

  if [ ! -f $OFILE  -o  ! -f $OFILE_REF ] ; then

#    echo Comparing $OFILE and $OFILE_REF":"
#    echo " "

#    for iono_var in Ne_max H_peak H_width ; do

#      var_cntl=$($NCDUMP -v$iono_var $OFILE_REF |grep "$iono_var = " |cut -d"=" -f2 |cut -d";" -f1 |sed -e 's/[eE]+*/\\*10\\^/')

#      var_test=$($NCDUMP -v$iono_var $OFILE     |grep "$iono_var = " |cut -d"=" -f2 |cut -d";" -f1 |sed -e 's/[eE]+*/\\*10\\^/')

#      if [ $var_cntl != $ropp_MDFV  -a  $var_cntl != $ropp_MDFV ] ; then

#        var_diff=$(echo "scale=10; $var_test - $var_cntl; " | bc -l)
#        var_diff_sq=$(echo "scale=10; $var_diff * $var_diff ; " | bc -l)

#        var_diff_flag=$(echo $var_diff_sq 1.0 | awk '{if ($1 > $2) print 1; else print 0}')

#        if [ "$var_diff_flag" -eq "1" ] ; then
#          echo "test $iono_var = " $var_test
#          echo "ref  $iono_var = " $var_cntl
#          echo "$iono_var in reference and test differs by more than 1.0 - test FAILED"
#          echo " "
#          ierr=$(($ierr + 1))
#        fi

#      fi

#    done

    echo "*** $OFILE or $OFILE_REF missing ***" >> $LOGFILE
    exit 1

  fi

# 7.2 Echo $LOGFILE if requested
# ------------------------------

  if [ -n ${ROPP_VERBOSE:-''} ] ; then
    cat $LOGFILE
  fi

# 7.3 Run ropp_1dvar_compare to compare results
# ---------------------------------------------

  ./ropp_1dvar_compare $OFILE $OFILE_REF $TEST $COMMENT


#-------------------------------------------------------------------------------
# 8. Finish
#-------------------------------------------------------------------------------

  echo "... examine $LOGFILE for details"
  echo " "


done  # loop over BGFILE


exit 0