#!/bin/bash
# Script reads daily concentrations and writes in SMAT-ready format
# Script relies on common linux commands and following programs
#
# proj -version: 5.0.1
# NCO  -version: 4.7.5
#
# Input: IOAPI file processed with hr2day
# Output: SMAT-ready CSV file
# User only needs to provide location of input and output file
# Options:
#
#    
#    Remove_Buffer: reads NX-2 and NY-2 columns and rows to remove
#                    buffer.
#
#    Include_Salt: write values for SALT. This was needed by MATS
#                   but SMAT needs to remove it. Recommended to not
#		   include salt for SMAT. Only matters for PM25 selection
#
#    format: C style printf format to increase or decrease output
#            decimal values for output
#
# this script expects to have the following species in the input 
# netcdf file. Names are hardwired
# PM25 (values in ug/m3): CRUSTAL NH4 SO4 EC NO3 OC PM25 CM
# O3   (values in ppb)  : O3

#####################################################################
# USER INPUTS


cases=( runS20_2026_0.9NOxred runS20_2026_0.9VOCred )
for scen in ${cases[@]}
do

echo "$scen"

infile=../1.hr2day/output/Denver_O3.${scen}.dailyavg.O3.4km.nc
outfile=/disk40/Denver_O3_2023/camx/postproc/3.SMAT/2.nc2SMAT/output/smat.denver23.input.o3.${scen}.csv
Remove_Buffer=No #Yes, set to "No" for CAMxv7.2 outputs since there is no buffer grid created (checked by ncview cell at ROW=COL=0 to find non-zero values)
SMAT_Type=O3 # Options: O3 or PM25
Include_Salt=No #Valid for PM25 only
format="%0.8f\n"

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

# FUNCTIONS

# Gets values from attributes
function ncattget { ncks --trd -M -m ${3} | grep -E -i "^${2} attribute [0-9]+: ${1}" | cut -f 11- -d ' ' | sort ; }

# Gets size of dimensions
function ncdmnsz { ncks --trd -m -M ${2} | grep -E -i "dimension .?: ${1}, size =" | cut -f 7 -d ' ' |uniq ; }

# Gets ID column for SMAT file
function create_ID {
ncols=${1}
nrows=${2}
ndays=${3}
xorig=${4}
yorig=${5}
xcell=${6}
ycell=${7}
echo $ycell
gdate0=${8}

echo "Ancillary variables..."
today=$gdate0
quote="\"\""

declare -a id
declare -a xdis
declare -a ydis
declare -a type
declare -a date

# Loop over columns and rows to get things that do not change in time
for (( r=1; r<=$nrows; r++ ))
do
	for (( c=1; c<=$ncols; c++ ))
	do

		tmp=$(( $c*1000+$r ))               #ID column
		id+=($tmp)


		x=$(echo "$xorig+$xcell/2+($c-1)*$xcell" | bc)  #XY for lat,lon, if xcell/ycell = 1333.333, e.g. float -> parse with |bc to do math
		y=$(echo "$yorig+$ycell/2+($r-1)*$ycell" | bc)
		xdis+=($x)
		ydis+=($y)

		type+=($quote)                      #Type column
		date+=($today)                      #Date column

	done #cols
done #rows

printf "%i\n" "${id[@]}">$TMP_DIR/tmp.id.1

printf "%.2f\n" "${xdis[@]}">$TMP_DIR/tmp.X.1
printf "%.2f\n" "${ydis[@]}">$TMP_DIR/tmp.Y.1
paste $TMP_DIR/tmp.X.1 $TMP_DIR/tmp.Y.1 >$TMP_DIR/tmp.XY.1
# Convert XY pairs to lat and lon columns
# This line of proj only works for LCC projections and assumes earth as a sphere.
proj -I -s +proj=lcc +a=6370997 +b=6370997 +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 -f '%.6f' <$TMP_DIR/tmp.XY.1 >$TMP_DIR/tmp.latlon.1
cut -f 1 $TMP_DIR/tmp.latlon.1 >$TMP_DIR/tmp.lat.1
cut -f 2 $TMP_DIR/tmp.latlon.1 >$TMP_DIR/tmp.lon.1

printf "%s\n" "${type[@]}">$TMP_DIR/tmp.type.1

printf "%i\n" "${date[@]}">$TMP_DIR/tmp.date.1

# Loop over days to update Date and increase  size of other ancillary columns
for (( day=1 ; day<=$ndays; day++))
do
	echo "Day: $today"
	yesterday=$today
	cat $TMP_DIR/tmp.id.1 >> $TMP_DIR/tmp.id
	cat $TMP_DIR/tmp.lat.1 >> $TMP_DIR/tmp.lat
	cat $TMP_DIR/tmp.lon.1 >> $TMP_DIR/tmp.lon

	cat $TMP_DIR/tmp.type.1 >> $TMP_DIR/tmp.type
	cat $TMP_DIR/tmp.date.1 >> $TMP_DIR/tmp.date

	today=$(date -d "$today +1 day" +%Y%m%d)
	sed -i "s/$yesterday/$today/g" $TMP_DIR/tmp.date.1

done # Days

}


function get_Data {
species=${1}
format=${2}
min=${3}
max=${4}
infile=${5}

ncks -H -s "$format" --trd -v $species -d COL,$min,$max -d ROW,$min,$max $infile |head -n -2 >>$TMP_DIR/tmp.$species

}

function create_Header {
subheader=${1}
echo "DAY" >$TMP_DIR/tmp.header
echo "$subheader" >>$TMP_DIR/tmp.header
}


#####################################################################
# MAIN

#create tmp folder
TMP_DIR=./tmp.SM.XY.$RANDOM
mkdir -p $TMP_DIR 

# Obtain values from file attributes
ncols=$(ncattget NCOLS global $infile)
nrows=$(ncattget NROWS global $infile)
ndays=$(ncdmnsz TSTEP $infile)

xorig=$(ncattget XORIG global $infile)
yorig=$(ncattget YORIG global $infile)

xcell=$(ncattget XCELL global $infile)
ycell=$(ncattget YCELL global $infile)

sdate=$(ncattget SDATE global $infile)
yyyy=${sdate:0:4}
jjj=${sdate:4}
gdate0=`date -d "$yyyy-01-01 + $(($jjj -1)) day" +%Y%m%d`

# Set indexes for NCO hyperslab (default is include buffer unless user removes it)
min=0
max=-1


# Remove buffer? Set variables
case $Remove_Buffer in
	Yes|yes|YES|TRUE|true|T|t)
	ncols=$(( $ncols-2 ))
	nrows=$(( $nrows-2 ))
	xorig=$(( $xorig + $xcell ))
	yorig=$(( $yorig + $ycell ))
	min=1
	max=-2
	;;
esac


# SMAT species are hardwired. 

# First CASE selects type of output: ozone or PM25(visibility)
case $SMAT_Type in
	Ozone|O3|ozone|o3)
	SMAT_species=(O3)
	header="_ID,_TYPE,LAT,LONG,DATE,O3"
	file_list="tmp.id tmp.type tmp.lat tmp.lon tmp.date tmp.O3"
	;;

        PM25|pm25|VIS|Vis|vis|Visibility|visibility)
	# Second case used to include or exclude SALT species
	case $Include_Salt in
		Yes|yes|YES|TRUE|true|T|t)
		SMAT_species=(CRUSTAL NH4 SO4 EC NO3 OC PM25 CM SALT)
		header="_ID,_TYPE,LAT,LONG,DATE,CRUSTAL,NH4,SO4,EC,NO3,OC,PM25,CM,SALT"
		file_list="tmp.id tmp.type tmp.lat tmp.lon tmp.date tmp.CRUSTAL tmp.NH4 tmp.SO4 tmp.EC tmp.NO3 tmp.OC tmp.PM25 tmp.CM tmp.SALT"
		;;

		No|no|NO|FALSE|false|F|f)
		SMAT_species=(CRUSTAL NH4 SO4 EC NO3 OC PM25 CM)
		header="_ID,_TYPE,LAT,LONG,DATE,CRUSTAL,NH4,SO4,EC,NO3,OC,PM25,CM"
		file_list="tmp.id tmp.type tmp.lat tmp.lon tmp.date tmp.CRUSTAL tmp.NH4 tmp.SO4 tmp.EC tmp.NO3 tmp.OC tmp.PM25 tmp.CM"
		;;

	esac
        ;;
esac




date

echo "Extracting concentrations..."
echo " "
for species in ${SMAT_species[@]}
do
	get_Data $species $format $min $max $infile
done

echo "Creating ancillary columns..."
echo " "
create_ID $ncols $nrows $ndays $xorig $yorig $xcell $ycell $gdate0 


echo "Create header and add to file..."
echo " "
create_Header $header
cd $TMP_DIR
paste -d "," $file_list > smat.outfile.data
cat tmp.header smat.outfile.data >smat.outfile
mv smat.outfile $outfile

echo "Succesful. Output: $outfile"
date

cd ../
rm -rf $TMP_DIR

done