#! /bin/bash

###------------------------------------------------------------------------------------------------###
###---                                                                                          ---###
###--- Realtime driver for DTOPS: Deterministic to Probabilistic Statistical Model (ECMWF ver.) ---###
###---                                                                                          ---###
###--- This script produces the realtime output for the DTOPS model.  It reads from the current ---###
###--- A-Deck file and calls Fortran routines to compute the DTOPS forecasts for each RI        ---###
###--- threshold.  The only required inputs are the TC ID (e.g., al07) and the full path to the ---###
###--- CARQ (.com) file. The user should also set the correct path to the realtime A-Deck files ---###
###--- on their system in the "Define paths, model codes, and coefficient files" section below. ---###
###---                                                                                          ---###
###--- Usage: $> ./get_dtops_realtime.bash al07 carq_comfile                                    ---###
###---                                                                                          ---###
###--- VERSION 1.3.1                                                                            ---###
###---                                                                                          ---###
###--- Edit history:                                                                            ---###
###---               Matt Onderlinde - May 2017 - wrote original code                           ---###
###---               Matt Onderlinde - Jun 2018 - updated to run on WCOSS                       ---###
###---               Stephanie Stevenson - Jul 2018 - updated to use merged a-deck file         ---###
###---               Stephanie Stevenson - Apr 2020 - updated for 2020 coef files, modified     ---###
###---  					      variable paths for WCOSS delivery         ---###
###---               Matt Onderlinde - Mar 2023 - updated from HWRF to HFSA/HFSB (HFSB not used ---###
###---                                            by DTOPS yet but plumbing is in place)        ---###
###---               Matt Onderlinde - Feb 2024 - updated for 2024 coefficients and text output ---###
###---                                            bug fix. Also added SHIPS Consensus+DTOPS     ---###
###---                                            consensus to text and E-deck output           ---###
###------------------------------------------------------------------------------------------------###

###--- Define paths, model codes, and coefficient files -------------------------------------------###

###--- Define directory for executables
EXECDIR="${HOMEnhcg}/exec/"

###--- Define A-Deck paths
adeck_path=${WORK_DIR}/storm-data/
#adeck_path=${ldir}/


###--- Define ATCF model codes
dshp_code=DSPE   ###--- ATCF model code for DSHP model
lgem_code=LGME   ###--- ATCF model code for LGEM model
hfai_code=HFAI   ###---	ATCF model code	for HFAI (HFSA interpolated) model
hfsa_code=HFSA   ###--- ATCF model code for HFSA (uninterpolated version)
hfbi_code=HFBI   ###--- ATCF model code for HFBI (HFSB interpolated) model
hfsb_code=HFSB   ###--- ATCF model code for HFSB (uninterpolated version)
avni_code=AVNI   ###---	ATCF model code	for AVNI (GFS interpolated) model
avno_code=AVNO   ###--- ATCF model cdoe for AVNO (GFS uninterpolated version)
emx_code=EMX     ###--- ATCF model code for EMX (ECMWF late) model

###--- Define DTOPS coefficient files
# Coefs and Data files ...
export DTOPS_COEF="$HOMEnhcg/parm/"
al_coefs_file=${DTOPS_COEF}/dtops_AL_coefs_2024.dat
ep_coefs_file=${DTOPS_COEF}/dtops_EP_coefs_2024.dat
al_emx_regress_file=${DTOPS_COEF}/dtops_AL_EMX_regress_2024.dat
ep_emx_regress_file=${DTOPS_COEF}/dtops_EP_EMX_regress_2024.dat

###--- Define output directory to which to send e-deck output files
outdir=${COMOUT}/dtops/

###--- Define output directory to which to send SHIPS-style text output files
text_outdir=${COMOUT}/dtops/

###--- Define directory where SHIPS output text files reside
ships_text_dir=${COMOUT}/ships/etext/

###--- Done defining paths, codes, and coefficient files ------------------------------------------###
###------------------------------------------------------------------------------------------------###
###------------------------------------------------------------------------------------------------###


###--- Get command line arguments -----------------------------------------------------------------###

tcid=${1}
carq_file=${2}
if [ ${#carq_file} -eq 0 ];then
	echo ""; echo "ERROR: Must pass in the TC ID and the full path to the CARQ (.com) file"
	echo "Usage: $> ./get_dtops_realtime.bash al07 /path/to/al902018.com"; echo ""
	exit
fi

###--- Done getting command line arguments --------------------------------------------------------###


###--- Make output directories if necessary -------------------------------------------------------###

mkdir -p -m 775 ${outdir}
mkdir -p -m 775 ${text_outdir}

###--- Done making output directories -------------------------------------------------------------###

###--- Define a function that will append final DTOPS output to the SHIPS text file
function append2ships {

	###--- Get inputs
	text_output=${1}
	text_outfile=${2}
	ships_text_file=${3}
	tcidyear=${4}
	eDeckOutFile=${5}
	vmax_init=${6}
	cyc=${7}

	###--- Report output to terminal
	echo "DTOPS output: "; echo "${text_output}"; echo ""

	###--- Construct SHIPS text file with DTOPS line added
	nlines_ships=`cat ${ships_text_file} | wc -l`
	consensus_line=`grep -n Consensus ${ships_text_file} | grep "[0-9]%" | head -1 | cut -f 1 -d ":"`
	n2tail=$(( ${nlines_ships} - ${consensus_line} ))
	head -${consensus_line} ${ships_text_file} > ${text_outfile}

	###--- If DTOPS couldnt run, fill with missing value
	if [ "${text_output}" == "MISSING" ];then
		echo "       DTOPS:   999.0%  999.0%  999.0%  999.0%  999.0%  999.0%  999.0%  999.0%" >> ${text_outfile}
		echo "       SDCON:   999.0%  999.0%  999.0%  999.0%  999.0%  999.0%  999.0%  999.0%" >> ${text_outfile}

	###--- Otherwise, write output to file
	else
		###--- Write DTOPS output to the file
		printf "%13s %8s%8s%8s%8s%8s%8s%8s%8s\n" ${text_output} >> ${text_outfile}

		###--- Compute SHIPS CONSENSUS + DTOPS consensus forecast (50% SHIPS-RII Consensus, 50% DTOPS); append to text/e-deck files
		fHours=(12 24 24 24 24 36 48 72)
		thresh=(20 25 30 35 40 45 55 65)
		read -a SCON <<< $( grep Consensus ${ships_text_file} | grep "[0-9]%" | head -1 | cut -f 2 -d ":" | sed -e 's/%//g' )
		read -a DTOP <<< $( echo "${text_output}" | cut -f 2 -d ":" | sed -e 's/%//g' )
		sdConText="       SDCON: "
		for t in $(seq 0 1 7)
		do
			SDCON[${t}]=$( echo "scale=1; (${SCON[t]} + ${DTOP[t]})/2.0" | bc )
			sdConText="${sdConText} ${SDCON[t]}%"
			int_plus_thresh=$(( ${vmax_init} + ${thresh[t]} ))
			tcidUpper=$( echo ${tcidyear:0:2} | tr '[:lower:]' '[:upper:]' )
			printf "%2s,%3s,%11s,%3s,%5s,%4s,%5s,%6s,%4s,%4s,%4s,%4s,%4s,%4s\n"     \
                               ${tcidUpper} ${tcidyear:2:2} ${cyc} RI SDCE ${fHours[t]} " " " " \
                               $( printf %.f ${SDCON[t]} ) ${thresh[t]} ${int_plus_thresh} " " 0 ${fHours[t]} >> ${eDeckOutFile}
		done
		printf "%13s %8s%8s%8s%8s%8s%8s%8s%8s\n" ${sdConText} >> ${text_outfile}
	fi

	###--- Add remaining SHIPS text file lines below
	tail -${n2tail} ${ships_text_file} >> ${text_outfile}

	###--- Cleanup
	tcid=${tcidyear:0:4}
	year=${tcidyear:4:4}
	rm ${tcid}${year}.avno ${tcid}${year}.avno.interp ${tcid}${year}.hfsa ${tcid}${year}.hfsa.interp \
           ${tcid}${year}.hfsb ${tcid}${year}.hfsb.interp ${tcid}${year}.hwrf ${tcid}${year}.hwrf.interp ${tcid}${year}.init_int ${tcid}${year}.int.inp 2> /dev/null
}

###--- Determine most recent synoptic cycle from which to compute DTOPS forecast ------------------###

cyc=`grep -i CARQ ${carq_file} | tail -1 | cut -f 3 -d "," | sed -e 's/ //g'`
if [ ${#cyc} -eq 0 ];then
	echo ""; echo "ERROR: Could not identify most recent model cycle... CARQ may be absent. Quitting here."; echo ""
	exit
fi
echo "Computing DTOPS RI percentages for ${cyc}"

###--- Done determining model cycle ---------------------------------------------------------------###


###--- Determine output e-deck file name ----------------------------------------------------------###

basin=`grep ${cyc} ${carq_file} | head -1 | cut -f 1 -d "," | sed -e 's/ //g' | tr '[:lower:]' '[:upper:]'`
snum=`echo ${tcid} | cut -c 3-4`
year=`echo ${cyc} | cut -c 1-4`
short_year=`echo ${year} | cut -c 3-4`
cyc_short=`echo ${cyc} | cut -c 3-10`
eDeckOutFile="${outdir}/EC_${cyc_short}${basin}${snum}${short_year}_dtops.edk"
text_outfile="${text_outdir}/EC_${cyc_short}${basin}${snum}${short_year}_dtops.txt"
ships_text_file="${ships_text_dir}/EC_${cyc_short}${basin}${snum}${short_year}_ships.txt"
tcidyear="${tcid}${year}"

###--- Done determining output file name ----------------------------------------------------------###


###--- Determine A-Deck file name -----------------------------------------------------------------###

adeck=`ls -d ${adeck_path}/${tcid}${year}/${tcid}${year}_interim_adeck.dat | tail -1`
if [ ${#adeck} -eq 0 ];then
	echo ""; echo "ERROR: Could not find A-Deck file. Quitting here."; echo ""
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
	exit
fi

###--- Done determining A-Deck file name ----------------------------------------------------------###


###--- Assemble an input A-deck from which to generate interpolated AVNI, HFAI, HFBI --------------###
###--- NOTE: This adeck will be used for all subsequent DTOPS generation             --------------###

###--- Get inputs
rm ${tcid}${year}.adeck 2> /dev/null
grep -i CARQ ${carq_file} | grep ", ${cyc}," > ${tcid}${year}.adeck
short_cycle=`echo ${cyc} | cut -c 3-10`
tcid_capitalized=`echo ${tcid} | tr '[:lower:]' '[:upper:]'`
ships_dat_file=${ships_text_dir}/EC_${short_cycle}${tcid_capitalized}${short_year}_ships.dat
aepoch=`date -d "${cyc:0:8} ${cyc:8:2}:00:00" +%s`
eepoch_minus6=$(( ${aepoch} - ( 3600 * 6 ) ))
prev_cyc=`date -d @${eepoch_minus6} "+%Y%m%d%H"`
grep ${prev_cyc} ${adeck}          | grep "${avno_code}" >> ${tcid}${year}.adeck
grep ${prev_cyc} ${adeck}          | grep "${hfsa_code}" >> ${tcid}${year}.adeck
grep ${prev_cyc} ${adeck}          | grep "${hfsb_code}" >> ${tcid}${year}.adeck
grep ${cyc}      ${ships_dat_file} | grep "${lgem_code}" >> ${tcid}${year}.adeck
grep ${cyc}      ${ships_dat_file} | grep "${dshp_code}" >> ${tcid}${year}.adeck
grep "${emx_code}" ${adeck} >> ${tcid}${year}.adeck
adeck=${tcid}${year}.adeck

###--- Grab initial intensity and forecast for HFSA, HFSB, AVNO
vmax_init=`grep ${cyc} ${adeck} | egrep "(CARQ,   0,|CARQ, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
if [ ${#vmax_init} -eq 0 ];then
        echo ""; echo "ERROR: Could not determine initial TC location (tried CARQ)... Quitting here..."; echo ""
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
        exit
fi
hfsa_check=`grep ${prev_cyc} ${adeck} | egrep "(${hfsa_code},   6,|${hfsa_code}, 006,)" | wc -l`
hfsb_check=`grep ${prev_cyc} ${adeck} | egrep "(${hfsb_code},   6,|${hfsb_code}, 006,)" | wc -l`
avno_check=`grep ${prev_cyc} ${adeck} | egrep "(${avno_code},   6,|${avno_code}, 006,)" | wc -l`
if [ ${hfsa_check} -eq 0 ] || [ ${avno_check} -eq 0 ];then
	echo ""; echo "ERROR: One or both of HFSA and AVNO is/are missing... Cant perform interpolation... Quiting here..."; echo ""
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
	exit
fi
rm ${tcid}${year}.init_int 2> /dev/null
echo "carq, ${vmax_init}" > ${tcid}${year}.init_int
hfsa_vinit=`grep ${prev_cyc} ${adeck} | egrep "(${hfsa_code},   6,|${hfsa_code}, 006,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
hfba_vinit=`grep ${prev_cyc} ${adeck} | egrep "(${hfsb_code},   6,|${hfsb_code}, 006,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
avno_vinit=`grep ${prev_cyc} ${adeck} | egrep "(${avno_code},   6,|${avno_code}, 006,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
echo "hfsa, ${hfsa_vinit}" >> ${tcid}${year}.init_int
echo "hfsb, ${hfsb_vinit}" >> ${tcid}${year}.init_int
echo "avno, ${avno_vinit}" >> ${tcid}${year}.init_int
rm ${tcid}${year}.hfsa ${tcid}${year}.hfsb ${tcid}${year}.avno 2> /dev/null
grep ${prev_cyc} ${adeck} | grep "${hfsa_code},...[0-9]," | cut -f 6,9 -d "," | sort | uniq | sort -n -t "," -k 1 > ${tcid}${year}.hfsa
grep ${prev_cyc} ${adeck} | grep "${hfsb_code},...[0-9]," | cut -f 6,9 -d "," | sort | uniq | sort -n -t "," -k 1 > ${tcid}${year}.hfsb
grep ${prev_cyc} ${adeck} | grep "${avno_code},...[0-9]," | cut -f 6,9 -d "," | sort | uniq | sort -n -t "," -k 1 > ${tcid}${year}.avno

###--- Call Fortran program to interpolate the HFSA/HFSB/AVNO forecasts
rm ${tcid}${year}.int.inp ${tcid}${year}.hfsa.interp ${tcid}${year}.hfsb.interp ${tcid}${year}.avno.interp 2> /dev/null
egrep -i "(CARQ|HFSA)" ${tcid}${year}.init_int > ${tcid}${year}.int.inp
${EXECDIR}/interpolate_adeck.x ${tcid}${year}.int.inp ${tcid}${year}.hfsa 0 120
egrep -i "(CARQ|HFSB)" ${tcid}${year}.init_int > ${tcid}${year}.int.inp
${EXECDIR}/interpolate_adeck.x ${tcid}${year}.int.inp ${tcid}${year}.hfsb 0 120
egrep -i "(CARQ|AVNO)" ${tcid}${year}.init_int > ${tcid}${year}.int.inp
${EXECDIR}/interpolate_adeck.x ${tcid}${year}.int.inp ${tcid}${year}.avno 6 120

###--- Append interpolated forecasts onto our new A-deck
tcid1=`echo ${tcid} | cut -c 1-2 | tr '[:lower:]' '[:upper:]'`
tcid2=`echo ${tcid} | cut -c 3-4 | tr '[:lower:]' '[:upper:]'`
spaces=" "
if   [ `echo "${vmax_init} < 10"  | bc` -eq 1 ];then
	spaces="   "
elif [ `echo "${vmax_init} < 100" | bc` -eq 1 ];then
	spaces="  "
fi
echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${hfai_code},   0, XXXX, XXXXX,${spaces}${vmax_init}, XXXX" >> ${adeck}
echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${hfbi_code},   0, XXXX, XXXXX,${spaces}${vmax_init}, XXXX" >> ${adeck}
echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${avni_code},   0, XXXX, XXXXX,${spaces}${vmax_init}, XXXX" >> ${adeck}
for n in 12 24 36 48 72
do
	###--- Define current forecast hour
	fhour=${n}

	###--- HFAI
	newint=""
	for t in $(seq 1 1 `cat ${tcid}${year}.hfsa.interp | wc -l`)
	do
		curfhr=`head -${t} ${tcid}${year}.hfsa.interp | tail -1 | cut -f 1 -d ","`
		if [ `echo "${curfhr} == ${n}" | bc` -eq 1 ];then
			newint=`head -${t} ${tcid}${year}.hfsa.interp | tail -1 | cut -f 2 -d "," | sed -e 's/ //g'`
			break
		fi
	done
	if [ ${#newint} -eq 0 ];then
		newint=0
	fi
	spaces=" "
	if   [ `echo "${newint} < 10"  | bc` -eq 1 ];then
		spaces="   "
	elif [ `echo "${newint} < 100" | bc` -eq 1 ];then
		spaces="  "
	fi
	echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${hfai_code},  ${fhour}, XXXX, XXXXX,${spaces}${newint}, XXXX" >> ${adeck}

	###--- HFBI
	newint=""
	for t in $(seq 1 1 `cat ${tcid}${year}.hfsb.interp | wc -l`)
	do
		curfhr=`head -${t} ${tcid}${year}.hfsb.interp | tail -1 | cut -f 1 -d ","`
		if [ `echo "${curfhr} == ${n}" | bc` -eq 1 ];then
			newint=`head -${t} ${tcid}${year}.hfsb.interp | tail -1 | cut -f 2 -d "," | sed -e 's/ //g'`
			break
		fi
	done
	if [ ${#newint} -eq 0 ];then
		newint=0
	fi
	spaces=" "
	if   [ `echo "${newint} < 10"  | bc` -eq 1 ];then
		spaces="   "
	elif [ `echo "${newint} < 100" | bc` -eq 1 ];then
		spaces="  "
	fi
	echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${hfbi_code},  ${fhour}, XXXX, XXXXX,${spaces}${newint}, XXXX" >> ${adeck}

	###--- AVNI
	newint=""
	for t in $(seq 1 1 `cat ${tcid}${year}.avno.interp | wc -l`)
	do
		curfhr=`head -${t} ${tcid}${year}.avno.interp | tail -1 | cut -f 1 -d ","`
		if [ `echo "${curfhr} == ${n}" | bc` -eq 1 ];then
			newint=`head -${t} ${tcid}${year}.avno.interp | tail -1 | cut -f 2 -d "," | sed -e 's/ //g'`
			break
		fi
	done
	if [ ${#newint} -eq 0 ];then
		newint=0
	fi
	spaces=" "
	if   [ `echo "${newint} < 10"  | bc` -eq 1 ];then
		spaces="   "
	elif [ `echo "${newint} < 100" | bc` -eq 1 ];then
		spaces="  "
	fi
	echo "${tcid1}, ${tcid2}, ${cyc}, XX, ${avni_code},  ${fhour}, XXXX, XXXXX,${spaces}${newint}, XXXX" >> ${adeck}
done

###--- Done assembling interpolated input A-deck --------------------------------------------------###


###--- Determine which EMX cycle to use and get initial minimum pressure --------------------------###

curhour=`echo ${cyc} | cut -c 9-10`
aepoch=`date -d "${cyc:0:8} ${cyc:8:2}:00:00" +%s`
eepoch_minus6=$(( ${aepoch} - ( 3600 * 6 ) ))
eepoch_minus12=$(( ${aepoch} - ( 3600 * 12 ) ))
ecyc6=`date  -d @${eepoch_minus6}  "+%Y%m%d%H"`
ecyc12=`date -d @${eepoch_minus12} "+%Y%m%d%H"`
emx_missing=0
if [ ${curhour} == "00" ] || [ ${curhour} == "12" ];then
	emx_pinit=`grep ${ecyc12} ${adeck} | grep "${emx_code},  12," | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
else
	p1=`   grep ${ecyc6}  ${adeck} | egrep "(${emx_code},   0,|${emx_code}, 000,)" | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
	p2=`   grep ${ecyc6}  ${adeck} | egrep "(${emx_code},  12,|${emx_code}, 012,)" | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
	if [ ${#p1} -gt 0 ] && [ ${#p2} -gt 0 ];then
		emx_pinit=`echo "scale=2; (${p1}+${p2})/2.0" | bc`
	fi
fi
if [ ${#emx_pinit} -eq 0 ];then
	echo ""; echo "Cannot find t=0 time in ECMWF A-deck... will replace with fitted value using other models..."; echo ""
	emx_missing=1
fi

###--- Done determining EMX model cycle and initial min P -----------------------------------------###


###--- Get initial intensities (VMAX) for DSHP, LGEM, HFAI, HFBI, AVNI ----------------------------###

dshp_vinit=`grep ${cyc} ${adeck} | egrep "(${dshp_code},   0,|${dshp_code}, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
lgem_vinit=`grep ${cyc} ${adeck} | egrep "(${lgem_code},   0,|${lgem_code}, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
hfai_vinit=`grep ${cyc} ${adeck} | egrep "(${hfai_code},   0,|${hfai_code}, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
hfbi_vinit=`grep ${cyc} ${adeck} | egrep "(${hfbi_code},   0,|${hfbi_code}, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
avni_vinit=`grep ${cyc} ${adeck} | egrep "(${avni_code},   0,|${avni_code}, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`

###--- Bail out if DSHP, LGEM, HFAI, HFBI, or AVNI is missing
core_missing=0
if [ ${#dshp_vinit} -eq 0 ];then
	core_missing=1
	echo "Could not find DSHP forecast for this cycle."
fi
if [ ${#lgem_vinit} -eq 0 ];then
	core_missing=1
	echo "Could not find LGEM forecast for this cycle."
fi
if [ ${#hfai_vinit} -eq 0 ];then
	core_missing=1
	echo "Could not find HFAI forecast for this cycle."
fi
#if [ ${#hfbi_vinit} -eq 0 ];then
#	core_missing=1
#	echo "Could not find HFBI forecast for this cycle."
#fi
if [ ${#avni_vinit} -eq 0 ];then
	core_missing=1
	echo "Could not find AVNI forecast for this cycle."
fi
if [ ${core_missing} -eq 1 ];then
	echo "Because at least one of DSHP, LGEM, HFAI, or AVNI was missing, cannot proceed."
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
	exit
fi

###--- Done getting initial intensities -----------------------------------------------------------###


###--- Get initial VMAX and initial TC latitude ---------------------------------------------------###

vmax_init=`grep ${cyc} ${adeck} | egrep "(OFCL,   0,|OFCL, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
tclat=`    grep ${cyc} ${adeck} | egrep "(OFCL,   0,|OFCL, 000,)" | head -1 | cut -f 7 -d "," | sed -e 's/ //g'`
tclonout=` grep ${cyc} ${adeck} | egrep "(OFCL,   0,|OFCL, 000,)" | head -1 | cut -f 8 -d "," | sed -e 's/ //g'`
if [ ${#vmax_init} -eq 0 ];then
	vmax_init=`grep ${cyc} ${adeck} | egrep "(CARQ,   0,|CARQ, 000,)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
fi
if [ ${#vmax_init} -eq 0 ];then
	echo ""; echo "ERROR: Could not determine initial VMAX (tried OFCL and CARQ)... Quitting here..."; echo ""
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
	exit
fi
if [ ${#tclat} -eq 0 ];then
	tclat=`   grep ${cyc} ${adeck} | egrep "(CARQ,   0,|CARQ, 000,)" | head -1 | cut -f 7 -d "," | sed -e 's/ //g'`
	tclonout=`grep ${cyc} ${adeck} | egrep "(CARQ,   0,|CARQ, 000,)" | head -1 | cut -f 8 -d "," | sed -e 's/ //g'`
fi
if [ ${#tclat} -eq 0 ];then
	echo ""; echo "ERROR: Could not determine initial TC location (tried OFCL and CARQ)... Quitting here..."; echo ""
	append2ships MISSING ${text_outfile} ${ships_text_file} ${tcidyear}
	exit
fi
tclatout=${tclat}
tclat=`echo ${tclat} | sed -e 's/[A-Z]//i'`
tclat=`echo "scale=2; ${tclat}/10.0" | bc`

###--- Done getting initial VMAX and TC latitude --------------------------------------------------###


###--- Determine ocean basin and corresponding coefficient files ----------------------------------###

if [ ${basin} == "AL" ];then
	cfile=${al_coefs_file}
	efile=${al_emx_regress_file}
else
	cfile=${ep_coefs_file}
	efile=${ep_emx_regress_file}
fi

###--- Done determining ocean basin ---------------------------------------------------------------###


###--- Loop through Rapid Intensification Thresholds and compute DTOPS percentages ----------------###

all_dt="12,24,24,24,24,36,48,72"
all_thresh="20,25,30,35,40,45,55,65"
text_output="       DTOPS: "
dshp_vfcst=${dshp_vinit}
lgem_vfcst=${lgem_vinit}
hfai_vfcst=${hfai_vinit}
hfbi_vfcst=${hfbi_vinit}
avni_vfcst=${avni_vinit}
emx_pfcst=${emx_pinit}
rm ${eDeckOutFile} 2> /dev/null
for i in $(seq 1 1 $(( `echo ${all_dt} | grep -o "," | wc -l` + 1 )))
do
	###--- Set current dt and threshold
	dt=`echo ${all_dt} | cut -f ${i} -d ","`
	thresh=`echo ${all_thresh} | cut -f ${i} -d ","`

	###--- Determine intensity (VMAX) for DSHP, LGEM, HFAI, HFBI, AVNI at the current forecast dt
	dshp_new_vfcst=`grep ${cyc} ${adeck} | egrep "(${dshp_code},  ${dt},|${dshp_code}, 0${dt},)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
	lgem_new_vfcst=`grep ${cyc} ${adeck} | egrep "(${lgem_code},  ${dt},|${lgem_code}, 0${dt},)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
	hfai_new_vfcst=`grep ${cyc} ${adeck} | egrep "(${hfai_code},  ${dt},|${hfai_code}, 0${dt},)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
	hfbi_new_vfcst=`grep ${cyc} ${adeck} | egrep "(${hfbi_code},  ${dt},|${hfbi_code}, 0${dt},)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
	avni_new_vfcst=`grep ${cyc} ${adeck} | egrep "(${avni_code},  ${dt},|${avni_code}, 0${dt},)" | head -1 | cut -f 9 -d "," | sed -e 's/ //g'`
	if [ ${#dshp_new_vfcst} -gt 0 ];then dshp_vfcst=${dshp_new_vfcst}; fi
	if [ ${#lgem_new_vfcst} -gt 0 ];then lgem_vfcst=${lgem_new_vfcst}; fi
	if [ ${#hfai_new_vfcst} -gt 0 ];then hfai_vfcst=${hfai_new_vfcst}; fi
	if [ ${#hfbi_new_vfcst} -gt 0 ];then hfbi_vfcst=${hfbi_new_vfcst}; fi
	if [ ${#avni_new_vfcst} -gt 0 ];then avni_vfcst=${avni_new_vfcst}; fi

	###--- Compute del-V for DSHP, LGEM, HFAI, HFBI, AVNI
	dshp_delv=`echo "scale=2; (${dshp_vfcst} - ${dshp_vinit})/1.0" | bc`
	lgem_delv=`echo "scale=2; (${lgem_vfcst} - ${lgem_vinit})/1.0" | bc`
	hfai_delv=`echo "scale=2; (${hfai_vfcst} - ${hfai_vinit})/1.0" | bc`
	hfbi_delv=`echo "scale=2; (${hfbi_vfcst} - ${hfbi_vinit})/1.0" | bc`
	avni_delv=`echo "scale=2; (${avni_vfcst} - ${avni_vinit})/1.0" | bc`

	###--- Determine del-P for EMX
	still_missing=0
	if [ ${emx_missing} -eq 0 ];then
		if [ ${curhour} == "00" ] || [ ${curhour} == "12" ];then
			emx_new_pfcst=`grep ${ecyc12} ${adeck} | egrep "(${emx_code},  $(( ${dt} + 12 )),|${emx_code}, 0$(( ${dt} + 12 )),)" | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
		else
			p1=`           grep ${ecyc6}  ${adeck} | egrep "(${emx_code},  ${dt},|${emx_code}, 0${dt},)"                         | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
			p2=`           grep ${ecyc6}  ${adeck} | egrep "(${emx_code},  $(( ${dt} + 12 )),|${emx_code}, 0$(( ${dt} + 12 )),)" | head -1 | cut -f 10 -d "," | sed -e 's/ //g'`
			if [ ${#p2} -gt 0 ];then
				emx_new_pfcst=`echo "scale=2; (${p1}+${p2})/2.0" | bc`
			else
				still_missing=1
			fi
		fi
		if [ ${still_missing} -eq 0 ];then
			if [ ${#emx_new_pfcst} -gt 0 ];then emx_pfcst=${emx_new_pfcst}; fi
			emx_delp=`echo "scale=2; (${emx_pfcst} - ${emx_pinit})/1.0" | bc`
		fi
	else
		emx_delp=`${EXECDIR}/fill_missing_EMX_realtime.x ${efile} ${dt} ${avni_delv} ${hfai_delv} ${dshp_delv} ${lgem_delv}`
	fi
	if [ ${still_missing} -eq 1 ];then
		echo "Had trouble determining EMX values for ${dt}... will fill using fitted value from other models."
		emx_delp=`${EXECDIR}/fill_missing_EMX_realtime.x ${efile} ${dt} ${avni_delv} ${hfai_delv} ${dshp_delv} ${lgem_delv}`
	fi

	###--- Compute DTOPS rapid intensification percentage
	rip=`${EXECDIR}/dtops_realtime.x ${cfile} ${dt} ${thresh} ${vmax_init} ${tclat} ${avni_delv} ${hfai_delv} ${dshp_delv} ${lgem_delv} ${emx_delp}`

	###--- Optional print-out of the result
	###  echo "dt, thresh, RI % = ${dt}, ${thresh}, ${rip}%"

	###--- Write output to E-Deck format
	int_plus_thresh=$(( ${vmax_init} + ${thresh} ))
	printf "%2s,%3s,%11s,%3s,%5s,%4s,%5s,%6s,%4s,%4s,%4s,%4s,%4s,%4s\n" \
               ${basin} ${snum} ${cyc} RI DTOP ${dt} " " " " ${rip} ${thresh} ${int_plus_thresh} " " 0 ${dt} >> ${eDeckOutFile}

	###--- Assemble probabilities output array
	if [ `echo "${rip} < 10" | bc` -eq 1 ];then
		sspaces="    "
	else
		sspaces="   "
	fi
	text_output="${text_output}${sspaces}${rip}.0%"

done

###--- Write DTOPS output text line (SHIPS-style) and append to SHIPS file
append2ships "${text_output}" ${text_outfile} ${ships_text_file} ${tcidyear} ${eDeckOutFile} ${vmax_init} ${cyc}