#!/bin/bash
export LC_ALL=C
ipipe=0
# Choose your solver here
# 0: direct skyline
# 1: parallel skyline
# 2: iterative sparse
solver=1


# Define log files
#-----------------
divarundir=$(pwd)
divacalclog=$divarundir'/output/divacalc.log'
divalog=$divarundir'/output/diva.log'

echo //////////////////////////////////////////////////  | tee $divacalclog
echo           Going to analyse a field                  | tee -a $divacalclog
echo //////////////////////////////////////////////////  | tee -a $divacalclog
#
divarunning=$(ps -eafl | grep -i 'diva\.a' | wc -l)
if [ "$divarunning" != "0" ]
then
echo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | tee -a $divacalclog
echo ????????diva already running?????? | tee -a $divacalclog
echo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | tee -a $divacalclog
ps -eafl | grep -i 'diva\.a'
fi
iqc=0
if [ "$1" == "qc" ] 
then
echo =============== >> $divacalclog
echo Quality control >> $divacalclog
echo =============== >> $divacalclog
iqc=1
fi

ibii=0
if [ "$1" == "Bii" ] 
then

echo =============== >> $divacalclog
echo Bii calculation >> $divacalclog
echo =============== >> $divacalclog
ibii=1
ipipe=0
fi


Filepar=./input/param.par
{
read linecomment
read lc
read linecomment
read icoordchange
read linecomment
read ispec
read linecomment
read ireg
read linecomment
read xori
read linecomment
read yori
read linecomment
read dx
read linecomment
read dy
read linecomment
read nx
read linecomment
read ny
read linecomment
read valex
read linecomment
read snr
read linecomment
read varbak
} < $Filepar

#######
# Stuff for errors
iflagcpme=0
iflagexerr=0
if [ "$ispec" -ge "100" ]
then
echo Clever poor mans error
iflagcpme=1
fi

if [ "$ispec" -le "-100" ]
then
if  [ "$ispec" -le "-110" ]
then
echo exerr with boundary effect
iflagexerr=2
ispec=17
else
echo exerr
iflagexerr=1
ispec=17
fi


fi

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





if [ "$ibii" == "1" ]
then
ispec=-1000
fi


if [ -f ./input/valatxy.coord ] 
then
cp ./input/valatxy.coord ./divawork/fort.79
fi
xi=1
cd ./meshgenwork
File23=./fort.23

cp  ../output/meshvisu/fort.10 .
cp  ../output/meshvisu/fort.22 .
cp  ../output/meshvisu/fort.23 .

{
read nnt1
read nnint
read nelt
} < $File23
 
cp fort.22 ../divawork/fort.11
cd ../divawork
# try to estimate memory allocation....
# 
NREA=0
NINT=0
NDATA=$(cat ../input/data.dat | wc -l)
echo 'Number of data points:' ${NDATA} >> $divacalclog
echo ' '
let NREA=$NREA+4*$NDATA
let NREA=$NREA+${nx}*${ny}
let NINT=$NINT+8*$NDATA
let NINT=$NINT+$nelt
let NINT=$NINT+41
let NDDLT=3*$nnt1+$nnint
let NDDLE=16
let NNT=$nnt1+$nnint
let NINT=$NINT+2*$NNT
let NINT=$NINT+$NDDLT
let NINT=$NINT+16*$nelt
let NINT=$NINT+2*$NDDLE
let NTERM=$NDDLT*$NDDLT
# instead of 90, square root of $NDDLT
let NDDLTB=$NDDLT
#echo Square root of $NDDLTB
#sqrt() { local x=$1 s=$1 os=0;while ((s!=os));\
#do os=$s;s=$(((s+x/s)/2)); done;echo $s; }
##sqrt() { local x=$1 s=$1 os=0; \
##while [ "$((os-s))" -ge "1" -o "$((s-os))" -ge "1" ]; \
##do os=$s;s=$(((s+x/s)/2)); done;echo $s; }
##xx=$(sqrt $NDDLTB)
xx=$(echo $NDDLTB | awk '{print int(sqrt($1))}')
#echo root is $xx $xxb
let NTERM=$NDDLT*$xx
let NTERM=$NTERM+$NTERM/2
let NREA=$NREA+$nelt
let NREA=$NREA+2*$NTERM
let NREA=$NREA+2*$NDDLT
let NREA=$NREA+$NDDLE*$NDDLE
let NREA=$NREA+$NDDLE*$NDDLE
let NREA=$NREA+$NDDLE
let NREA=$NREA+38*$nelt
let NREA=$NREA+2*$nnt1
# For optimisation of element location, rought estimate:
let NINT=$NINT+$nelt*100
#
itcs=0
if [ -f ../input/constraint.dat ]
then
let itcs=$itcs+1

echo Advection requested >> $divacalclog

cat ../input/constraint.dat > ./fort.49
cp ../input/UVinfo.dat ./fort.93
cp ../input/Uvel.dat ./fort.91
cp ../input/Vvel.dat ./fort.92

#JMB2012 added support for source terms in advection constraint
# as many sources as whished, format: x y Q
if [ -f ../input/sources.dat ]
then
cp ../input/sources.dat ./fort.24
fi


let NREA=$NREA+2*$nnt1
{
read aaaaa
read bbbbb
read ccccc
read ddddd
read nxuuu
read nyuuu
} < ../input/UVinfo.dat
let NREA=$NREA+$nxuuu*$nyuuu
let NREA=$NREA+2*$nnt1
let NREA=$NREA+1
fi
if [ -f ../input/RL.dat ]
then
let itcs=$itcs+2
echo Variable L requested >> $divacalclog
cp ../input/RL.dat ./fort.94
cp ../input/RLinfo.dat ./fort.95
let NREA=$NREA+$nnt1
{
read aaaaa
read bbbbb
read ccccc
read ddddd
read nxlll
read nylll
} < ../input/RLinfo.dat
let NREA=$NREA+$nxlll*$nylll
let NREA=$NREA+$nnt1
let NREA=$NREA+1
fi

#creation of file fort.13
#
# 'griddef.a' replaced by awk command 
#------------------------------------
#echo $xori > ./fort.13
#echo $yori >> ./fort.13
#echo $dx >> ./fort.13
#echo $dy >> ./fort.13
#echo $nx >> ./fort.13
#echo $ny >> ./fort.13
#echo $valex >> ./fort.13
#../../bin/griddef.a >> diva.log

echo $(echo $xori $dx | awk '{print ($1-$2)}') >  ./fort.13
echo $(echo $yori $dy | awk '{print ($1-$2)}') >> ./fort.13
echo $dx >> ./fort.13
echo $dy >> ./fort.13
echo $nx >> ./fort.13
echo $ny >> ./fort.13
echo $valex >> ./fort.13


echo $(echo $lc | awk '{print (1/($1*$1*$1*$1))}') > fort.12
echo $(echo $lc $xi | awk '{print (2*$2/($1*$1))}') >> fort.12

if [ "$ireg" == "3" ]
then
aaaref=$(head -1 fort.12)
echo 0 > newalpha
# Test, not really zero, but small

##echo  1.E-6 $aaaref | ../../bin/multiply.a > newalpha
echo $(echo $aaaref | awk '{print ($1*1.E-6)}') >  newalpha

tail -1 fort.12 >> newalpha
mv newalpha fort.12
let ireg=1
fi
echo $varbak > ./fort.15
cd ..
if [ "$iqc" != "1" ]
then
if [ "$ispec" -lt "0" ] 
then
if [ "$varbak" != "0" ]
then
if [ "$1" != "-pfsum" ]
then
#dvkern
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog
echo %  Going to call another diva execution    % >> $divacalclog
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog


echo divacalc: ////////////////////////////////////////// >> $divalog
echo divacalc:         running dvcovariance >> $divalog
echo divacalc: ////////////////////////////////////////// >> $divalog
echo ' '  >> $divalog

. dvcovariance
    cat $divacalclog ./output/dvcovariance.log > bidon
    mv bidon $divacalclog

if [ "$ibii" == "1" ]
then
cd divawork
#pwd
../../bin/bintoasc.a | awk '{print $1,$2,$3/(1-$3)}' >../output/Bii.dat
cd ..
echo Finished bii  >> $divacalclog
exit
fi

if [ "$ipipe" == "1" ]
then
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog
echo %  Have launched another diva execution    % >> $divacalclog
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog
else
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog
echo %  Have finished another diva execution    % >> $divacalclog
echo %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% >> $divacalclog
fi
fi
fi
fi
fi
if [ -f ./input/data.dat ] 
then 
cp ./input/data.dat ./divawork/fort.44
else
echo Need to provide data.dat in ./input !
##ls -l ./input
##pwd
fi
cd ./divawork
echo coord 1 > ./fort.10
echo $icoordchange >> ./fort.10
echo mathpr 1 >> ./fort.10
echo 2 >> ./fort.10
echo 0 >> ./fort.10
echo 2 >> ./fort.10
echo topolo 1 >> ./fort.10
echo $nnt1 >> ./fort.10
echo $nnint >> ./fort.10
echo $nelt >> ./fort.10
echo datapr 1 >> ./fort.10
echo $ireg >> ./fort.10
#JMB
if [ "$itcs" != "0" ]
then
echo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! >> $divacalclog
echo  Module constraint activated >> $divacalclog
echo $itcs
echo constr 1 >> ./fort.10
echo $itcs >> ./fort.10
fi
echo solver 1 >> ./fort.10
echo $solver >> ./fort.10
if [ "$solver" == "2" ]
then
echo 0 0 0 >> ./fort.10
fi
echo stores 1 >> ./fort.10
echo 3 >> ./fort.10
echo VARBAK: $varbak >> $divacalclog
if [ "$iqc" != "1" ]
then
if [ "$varbak" != "0" ]
then
if [ "$ispec" != "0" ]
then
	echo ' ' >> $divacalclog
echo Errors will be calculated >> $divacalclog
echo esterr 1 >> ./fort.10
echo $ispec >> ./fort.10
fi
fi
#echo $ >> ./fort.10
echo gcvfac 1 >> ./fort.10
echo 5 >> ./fort.10
echo stopex >> ./fort.10
else
#QC
echo dataqc 1 >> ./fort.10
echo stopex >>./fort.10
fi
#creation of file fort.20
nbccol=$(head -n 1 ./fort.44 | wc -w)
#echo $lc $snr $nbccol $xi | ../../bin/calcmu.a >> diva.log

# New version
multconst=$(echo $lc $snr $xi | ../../bin/calcmuc.a)
multvar=$(echo $multconst | awk '{print 1*$1}')

#echo multiply weight $multvar $multconst so what
awk -v multconst=$multvar '{ if (NF>2) {if ($4 > "") {$4=$4*multconst} else {$4=multconst}}; {print $1,$2,$3,$4}}' fort.44 > fort.20

##?? if [ $? -ne 0 ];then
##?? echo ' '
##?? echo --------------------------------------------
##?? echo A problem was encountered during execution !
##?? echo       divacalc:  calmu.a
##?? echo Check execution track
##?? echo --------------------------------------------

NDATAB=$(cat ./fort.20 | wc -l)
if [ "$NDATAB" -lt "$NDATA" ]
then
echo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | tee -a $divacalclog
echo ------------------------------------------- | tee -a $divacalclog
echo ??????????????????????????????????????????? | tee -a $divacalclog
echo ?????   Corrupted data file    ???????????? | tee -a $divacalclog
echo From $NDATA lines in data.dat 		 | tee -a $divacalclog
echo Only $NDATAB have been feed into DIVA 	 | tee -a $divacalclog
echo ??????????????????????????????????????????? | tee -a $divacalclog
echo ------------------------------------------- | tee -a $divacalclog
echo !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! | tee -a $divacalclog
fi
if [ "$1" == "-pfsum" ]
then
echo Summing covariances >> $divacalclog
cd ..

echo divacalc: ////////////////////////////////////////// >> $divalog
echo divacalc:         running dvcovariancesum >> $divalog
echo divacalc: ////////////////////////////////////////// >> $divalog
echo ' '  >> $divalog
. ./dvcovariancesum
    cat $divacalclog ./output/dvcovariancesum.log > bidon
    mv bidon $divacalclog
exit
fi

echo divacalc: ////////////////////////////////////////// >> $divacalclog
echo divacalc:'        Diva calculation (diva.a)'  >> $divacalclog
echo divacalc: ////////////////////////////////////////// >> $divacalclog
echo ' '  >> $divacalclog

##echo 'REMOVE FORT.35 and FORT.43'
rm -f fort.43
rm -f fort.35
#execute diva
echo $NINT $NREA | ../../bin/diva.a >> $divacalclog

if [ $? -ne 0 ];then
echo ' ' | tee -a $divacalclog $divalog
echo ------------------------------------------- | tee -a $divacalclog $divalog
echo A problem was encountered during execution !| tee -a $divacalclog $divalog
echo       divacalc:  diva.a			 | tee -a $divacalclog $divalog
echo 'Check execution track in diva.log'	 | tee -a $divacalclog $divalog
echo --------------------------------------------| tee -a $divacalclog $divalog
fi

dataonmesh=0
##ctroupin Oct 2012
# Check if fort.84 was generated AND if not empty...
if [ ! -s ./fort.84 ];
then
        echo ' ' | tee -a $divacalclog
        echo '//////////////////////////////////////////////////' | tee -a $divacalclog $divalog
        echo 'Problem in diva.a : analysis could not be finished' | tee -a $divacalclog $divalog
        echo ' '						  | tee -a $divacalclog $divalog
        echo 'Check diva.log'					  | tee -a $divacalclog $divalog
        echo '//////////////////////////////////////////////////' | tee -a $divacalclog $divalog
fi
##---




#if [ -f ./fort.35  ]; then
#dataonmesh=$(cat ./fort.35 | awk '{print $1-$2}')
#echo Data are outside the mesh, will produce empty files
#fi
#if [ "$dataonmesh" -lt "1" ]; then
#cd ..
#dvnil
#cd ./divawork
#else
#if [ -f ./fort.43 ] 
#then
#echo No data, will produce empty files
#echo $nx $ny $valex | ../../bin/ghervalex.a
#cp fort.84 fort.87
#fi
#fi

# New version 11 Jan 2012
#------------------------

if [ -f ./fort.35  ]; then
	dataonmesh=$(cat ./fort.35 | awk '{print $1-$2}')
fi
if [ "$dataonmesh" -lt "1" ]; then
	cd ..
	echo Data are outside the mesh, will produce empty files
	dvnil
	cd ./divawork
fi

rm -f fort.59
rm -f fort.61
if [ "$iqc" == "1" ] 
then
#
#output of results to user
cp fort.76 ../output/expectederroratdatapoint.anl
cp fort.71 ../output/fieldatdatapoint.anl


echo divacalc: ////////////////////////////////////////// | tee -a $divalog $divacalclog
echo divacalc:       lookforoutliers  			  | tee -a $divalog $divacalclog 
echo divacalc: ////////////////////////////////////////// | tee -a $divalog $divacalclog
echo ' '  						  | tee -a $divalog $divacalclog

echo $valex | ../../bin/lookforoutliers.a >> $divalog
if [ $? -ne 0 ];then
echo ' '
echo --------------------------------------------
echo A problem was encountered during execution !
echo       divacalc:   lookforoutliers.a
echo 'Check execution track in diva.log'
echo --------------------------------------------
echo ' ' >> $divalog
echo -------------------------------------------- >> $divalog
echo A problem was encountered during execution ! >> $divalog
echo       divacalc:    lookforoutliers.a  >> $divalog
echo 'Check execution track in diva.log' >> $divalog
echo -------------------------------------------- >> $divalog
fi
echo ' '

echo ' ' >> $divacalclog
echo Copying outlier lists >> $divacalclog
echo ' '
if [ -f fort.66 ]
then
mv -f fort.66 ../output/outliers.dat
mv -f fort.67 ../output/outliers.normalized.dat
else
echo ' '                                          | tee -a $divacalclog
echo -------------------------------------------- | tee -a $divacalclog
echo A problem was encountered during execution ! | tee -a $divacalclog
echo 'Check execution track in diva.log'	  | tee -a $divacalclog
echo -------------------------------------------- | tee -a $divacalclog
fi
else
#output as netcdf not yet done
echo ' ' >> $divacalclog
echo Output of results copied in ./output/ >> $divacalclog
echo ' ' >> $divacalclog
cp fort.84 ../output/fieldgher.anl
if [ -f fort.82 ]
then
cp fort.82 ../output/valatxyascii.anl
fi
cp fort.83 ../output/fieldascii.anl
if [ -f fort.87 ]
then
cp fort.87 ../output/errorfieldgher.anl
fi
if [ -f fort.73 ]
then
cp fort.73 ../output/erroratxyascii.anl
fi
if [ -f fort.86 ]
then
cp fort.86 ../output/errorfieldascii.anl
fi
cp fort.71 ../output/fieldatdatapoint.anl
if [ -f fort.72 ] 
then
cp fort.72 ../output/erroratdatapoint.anl
fi
if [ -f fort.77 ]
then
cp fort.77 ../output/gcvval.dat
fi

##echo ' '
##echo Creation of file GridInfo.dat
##echo ' '

echo $xori > ../output/ghertonetcdf/GridInfo.dat
echo $yori >> ../output/ghertonetcdf/GridInfo.dat
echo $dx >> ../output/ghertonetcdf/GridInfo.dat
echo $dy >> ../output/ghertonetcdf/GridInfo.dat
echo $nx >> ../output/ghertonetcdf/GridInfo.dat
echo $ny >> ../output/ghertonetcdf/GridInfo.dat
echo $valex >> ../output/ghertonetcdf/GridInfo.dat


#Iceedge case
if [ -f ../input/iceedge.dat ]
then
echo Ice edge to be used as mask
cd ..
dvicemask
cp ./output/fieldgher.anl ./divawork/fort.84
cd divawork
echo finished ice edge case
fi
cp fort.84 ../output/ghertonetcdf
if [ -f fort.87 ]
then 
cp  fort.87 ../output/ghertonetcdf
fi
iecal=0
if [ "$ispec" == "1" ]
then
iecal=1
fi
if [ "$ispec" == "3" ]
then
iecal=1
fi
if [ "$ispec" == "5" ]
then
iecal=1
fi
if [ "$ispec" == "7" ]
then
iecal=1
fi
if [ "$ispec" == "11" ]
then
iecal=1
fi
if [ "$ispec" == "13" ]
then
iecal=1
fi
if [ "$ispec" == "15" ]
then
iecal=1
fi
if [ "$ispec" == "17" ]
then
iecal=1
fi
if [ "$ispec" == "-1" ]
then
iecal=1
fi
if [ "$ispec" == "-3" ]
then
iecal=1
fi
if [ "$ispec" == "-5" ]
then
iecal=1
fi
if [ "$ispec" == "-7" ]
then
iecal=1
fi
if [ "$ispec" == "-11" ]
then
iecal=1
fi
if [ "$ispec" == "-13" ]
then
iecal=1
fi
if [ "$ispec" == "-15" ]
then
iecal=1
fi
if [ "$ispec" == "-17" ]
then
iecal=1
fi
# JMB ADAPTED FOR DIFFERENT CASES IF ISPEC
if [ "$varbak" != "0" ]
then
if [ "$iecal" == "1" ] 
then
echo Creating netcdf file for field and associated error >> $divacalclog
echo in directory ./output/ghertonetcdf/ >> $divacalclog
cd ../output/ghertonetcdf

echo divacalc: ////////////////////////////////////////// >> $divalog
echo divacalc:       netcdfoutput   >> $divalog
echo divacalc: ////////////////////////////////////////// >> $divalog
echo ' '  >> $divalog

../../../bin/netcdfoutput.a >> $divalog
if [ $? -ne 0 ];then
echo ' ' | tee -a $divacalclog $divalog
echo -------------------------------------------- | tee -a $divacalclog $divalog
echo A problem was encountered during execution ! | tee -a $divacalclog $divalog
echo       divacalc:   netcdfoutput.a		  | tee -a $divacalclog $divalog
echo 'Check execution track in diva.log'	  | tee -a $divacalclog $divalog
echo -------------------------------------------- | tee -a $divacalclog $divalog
fi
fi
fi
if [ "$varbak" == "0"  -o  "$iecal" == "0" ] 
then
echo Creating netcdf file only for field  >> $divacalclog
echo since Varbak and ispec are $varbak $ispec >> $divacalclog

cd ../output/ghertonetcdf
../../../bin/netcdfoutputfield.a >> $divalog
if [ $? -ne 0 ];then
echo ' ' | tee -a $divacalclog $divalog
echo -------------------------------------------- | tee -a $divacalclog $divalog
echo A problem was encountered during execution ! | tee -a $divacalclog $divalog
echo       divacalc:   netcdfoutputfield.a        | tee -a $divacalclog $divalog
echo 'Check execution track in diva.log'          | tee -a $divacalclog $divalog
echo -------------------------------------------- | tee -a $divacalclog $divalog
echo ' '                                          | tee -a $divacalclog $divalog
fi
###rm -f ../error*.anl
fi
cd ../..

echo ' ' | tee -a $divacalclog
echo -------------------- | tee -a $divacalclog
echo Analysis is finished | tee -a $divacalclog
echo -------------------- | tee -a $divacalclog


if [ "$iflagcpme" == "1" ]
then
divacpme 
fi

if [ "$iflagexerr" == "1" ]
then
divaexerr
fi


if [ "$iflagexerr" == "2" ]
then
divaexerr full
fi






fi
if [ "$ispec" -lt "0" ] 
then
if [ "$varbak" != "0" ]
then
rm -f ./divawork/sub/fort.*
fi
fi


