#!/bin/sh
set +xa
echo "---------------------------------------------------------------"
echo "Script - SSTOIQD PREP processing (NCDC quarter-degree SST analysis)"
echo "---------------------------------------------------------------"
echo "History: October 2010 - Implement SSTOIQD Analysis"
echo "---------------------------------------------------------------"
set -xa

###########################################################
# Processing/Flow of script
# 1) Copy this job's input files to the working directory
# 2) For each date used in satellite bias correction (and analysis):
# 2a)   - Generate buoy and ship superobs 
# 2b)   - Apply bias correction to ship superobs 
# 2c)   - Generate avhrr superobs (for each avhrr platform used/available)
# 2d)   - Generate amsr superobs (if amsr data used/available)
# 3) Save $PDYm1 buoy superobs and corrected ship superobs to $COMOUT.
# 4) Generate ice superobs:
#     - Compute median concentration over $ndays_ice days
#     - Generate simulated SSTs from ice concentration values.
# 4) Save ice superobs to $COMOUT.
# 5) For each satellite used and for day vs night:
# 5a) Compute satellite bias correction weights
# 5b) Apply bias correction to satellite superobs for $PDYm1
# 5c) Save $PDYm1 corrected satellite superobs to $COMOUT.
# 7) Protect restricted ship files sent to $COMOUT
###########################################################

#############################################################
# Ensure that this job is pointed to the production temporary
# working directory, which is where all of this job's input 
# and output files will be saved during execution
#############################################################
cd $DATA

msg="HAS BEGUN on `hostname`"
postmsg "$jlogfile" "$msg"

AVHRR_name=${AVHRR_name:-sstnvh}           # bufr dumplist name for AVHRR data used
AMSR_name=${AMSR_name:-amsre}              # amsr data used
AVHRR_sats=${AVHRR_sats:-"noaa19 metop2"}  # avhrr satellites used
ndays_satbc=${ndays_satbc:-7}              # number of days used in computing satellite bias correction
ndays_ice=${ndays_ice:-7}                  # number of days used to compute ice field

# Executables
INSITUX=${INSITUX:-$EXECsstoiqd/sstoiqd_insitu}        # executable to generate ship/buoy superobs
AVHRRX=${AVHRRX:-$EXECsstoiqd/sstoiqd_avhrr}           # executable to generate AVHRR superobs
AMSRX=${AMSRX:-$EXECsstoiqd/sstoiqd_amsr}              # executable to generate AMSR superobs
SHIPBCX=${SHIPBCX:-$EXECsstoiqd/sstoiqd_shipbias4}     # executable to bias correct ship superobs
ICEGRDX=${ICEGRDX:-$EXECsstoiqd/sstoiqd_ice2_1440x720} # executable to regrid ice data 
ICE2SSTX=${ICE2SSTX:-$EXECsstoiqd/sstoiqd_ice2sst_med} # executable to generated simulated SST superobs from ice data.
SATBWTX=${SATBWTX:-$EXECsstoiqd/sstoiqd_eotbiaswt4}    # executable to generate sat bias weights
SATBCX=${SATBCX:-$EXECsstoiqd/sstoiqd_eotbiascor4}     # executable to apply bias correction to sat superobs

# Fixed files
MASK=${MASK:-$FIXsstoiqd/sstoiqd_quarter-mask-extend}             # quarter-degree mask
STDEV=${STDEV:-$FIXsstoiqd/sstoiqd_stdev1d-coads3-fill}           # obs standard deviation field
ICEMASK=${ICEMASK:-$FIXsstoiqd/sstoiqd_ice_flags_mask.dat}        # ice mask
ICECOEF=${ICECOEF:-$FIXsstoiqd/sstoiqd_gsfc-fit-coef-fill-final}  # ice to sst fit coefficients
FMODE=${FMODE:-$FIXsstoiqd/sstoiqd_eot6.damp-zero.ev130.ano.dat}  # EOT modes file (filled)
MASK2D=${MASK2D:-$FIXsstoiqd/sstoiqd_lstags.twodeg.dat}           # EOT Land/SEA FILE
CLIM2D=${CLIM2D:-$FIXsstoiqd/sstoiqd_clim.71.00.gdat}             # EOT Climate FILE
VARMD=${VARMD:-$FIXsstoiqd/sstoiqd_var-mode}                      # Variance of each modes

#############################################################
# cp input files to working directory
#############################################################

SATBC_DATElne_minus=`sh $USHutil/finddate.sh $PDY s-$ndays_satbc`
SATBC_DATElne_plus=`sh $USHutil/finddate.sh $PDY s+$ndays_satbc`

set +x
echo " "
echo "#############################################################"
echo "Dates considered for sat bias correction: "
echo "  $SATBC_DATElne_minus $PDY $SATBC_DATElne_plus"
echo "Any insitu or satellite files in $COMIN with these dates will be copied to working directory."
echo "Source codes will decide which dates to use."
echo "#############################################################"
echo " "
set -x

for CDATE in $SATBC_DATElne_minus $PDY $SATBC_DATElne_plus
do
 # bufr data
   for CTYPE in ships dbuoy mbuoy  ${AVHRR_name} ${AMSR_name}
   do
     file=$COMIN/sstoiqd.$CTYPE.$CDATE.bufr_d
     [ -s $file ] && cp $file .
   done
 # gridded amsr fields, if available...
   file=$COMIN/sstoiqd.${AMSR_name}grid.$CDATE.nc
   [ -s $file ] && cp $file .
 # guess file for qc of the data
   file=$COMIN/sstoiqd.sst.$CDATE.grb
   [ -s $file ] && cp $file .
done


ICE_DATElne_minus=`sh $USHutil/finddate.sh $PDY s-$ndays_ice`

set +x
echo " "
echo "#############################################################"
echo "Dates considered for computing median ice field: "
echo "  $ICE_DATElne_minus "
echo "Any ice files in $COMIN with these dates will be copied to working directory and processed."
echo "#############################################################"
echo " "
set -x

for CDATE in $ICE_DATElne_minus 
do
   file=$COMIN/sstoiqd.ice.$CDATE.grb
   [ -s $file ] && cp $file .
done


# settings for data coverage warnings
cov_buoy=0.1
cov_ship=0.05
cov_buoyship=0.0
cov_sat=2.5
cov_sat2=0.0
cov_amsr=15.0
cov_amsr2=1.0
cov_ice=6.0
cov_ice2=0.0

for CDATE in $SATBC_DATElne_minus $PDY $SATBC_DATElne_plus
do

iyrmid=`echo $CDATE | cut -c 1-4`
immid=`echo $CDATE | cut -c 5-6`
idmid=`echo $CDATE | cut -c 7-8`

# Process In Situ data

set +x
echo " "
echo "############# generate buoy ship superobs for $CDATE. ########################"
set -x

# define the input file 
#  later, process subtypes separately.  requires significant changes to NCDC source code.
#  for now, cat ship, dbuoy, mbuoy into one file.

flist=combfr.$CDATE.lst
[ -f $flist ] && rm $flist
for CTYPE in ships dbuoy mbuoy
do 
  [ -s sstoiqd.$CTYPE.$CDATE.bufr_d ] && echo sstoiqd.$CTYPE.$CDATE.bufr_d >> $flist
done
if [ -s $flist ]; then
  pgm=bufr_combfr
  . prep_step

  export XLFUNIT_50="insitu.$CDATE.bufr_d"
  $EXECbufr/bufr_combfr < $flist >> $logfile
  msg="Done combining insitu data for $CDATE"
  postmsg "$jlogfile" "$msg"
fi

[ -s insitu.$CDATE.bufr_d ] || continue

pgm=sstoiqd_insitu
. prep_step

export XLFUNIT_21="insitu.$CDATE.bufr_d"
export XLFUNIT_31="sstoiqd.sst.$CDATE.grb"
export XLFUNIT_30="$MASK"
export XLFUNIT_34="$STDEV"
export XLFUNIT_51="buoy.obs.$CDATE.txt"
export XLFUNIT_52="ship.obs.$CDATE.txt"
export XLFUNIT_61="buoy.grid.$CDATE.ieee"
export XLFUNIT_62="ship.grid.$CDATE.ieee"
export XLFUNIT_71="buoy.sobs.$CDATE.ieee"
export XLFUNIT_72="ship.sobs.$CDATE.ieee"

export XLFUNIT_75="ship.txt.$CDATE"
export XLFUNIT_76="moored.txt.$CDATE"
export XLFUNIT_77="drift.txt.$CDATE"

cat <<oiEOF > buoyship.parm
$iyrmid $immid $idmid
$cov_buoy
$cov_ship
$cov_buoyship
oiEOF

startmsg
$INSITUX <buoyship.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating insitu superobs for $CDATE"
postmsg "$jlogfile" "$msg"

set +x
echo " "
echo "############# apply bias correction to ship superobs for $CDATE. ########################"
set -x

pgm=sstoiqd_shipbias4
. prep_step

export XLFUNIT_12="ship.sobs.$CDATE.ieee"
export XLFUNIT_52="ship.cobs.$CDATE.ieee"

globc=0.14                              # global ship bias adjustment
cat <<inxEOF> shipbias4.inp
$globc
inxEOF

$SHIPBCX <shipbias4.inp >> $pgmout 2>errfile
msg="Done applying bias correction to ship superobs for $CDATE"
postmsg "$jlogfile" "$msg"


# Process AVHRR Data

for SATname in $AVHRR_sats; do

set +x
echo " "
echo "############# generate $SAT AVHRR superobs for $CDATE. ########################"
set -x

case $SATname in
  noaa17) isource=6; isaid=208;;
  noaa18) isource=7; isaid=209;;
  noaa19) isource=8; [ $AVHRR_name = sstns ] && isource=9; isaid=223;;
  metop2) isource=12; [ $AVHRR_name = sstns ] && isource=8; isaid=004;;
esac

pgm=sstoiqd_avhrr
. prep_step

export XLFUNIT_21="sstoiqd.$AVHRR_name.$CDATE.bufr_d"
export XLFUNIT_31="sstoiqd.sst.$CDATE.grb"
export XLFUNIT_30="$MASK"
export XLFUNIT_34="$STDEV"
export XLFUNIT_51="${SATname}day.obs.$CDATE.txt"
export XLFUNIT_52="${SATname}nte.obs.$CDATE.txt"
export XLFUNIT_61="${SATname}day.grid.$CDATE.ieee"
export XLFUNIT_62="${SATname}nte.grid.$CDATE.ieee"
export XLFUNIT_71="${SATname}day.sobs.$CDATE.ieee"
export XLFUNIT_72="${SATname}nte.sobs.$CDATE.ieee"

cat <<oiEOF > satellite.parm
$iyrmid $immid $idmid
$isaid
$cov_sat
$cov_sat2
oiEOF

startmsg
$AVHRRX <satellite.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating $SATname superobs for $CDATE"
postmsg "$jlogfile" "$msg"

done  # end of avhrr satellites loop


# Process AMSR Data

if [ -s sstoiqd.${AMSR_name}grid.$CDATE.nc ]; then

set +x
echo " "
echo "############# generate AMSR superobs for $CDATE. ########################"
set -x

pgm=sstoiqd_amsr
. prep_step

set +x
echo No XLFUNIT for sstoiqd.${AMSR_name}grid.$CDATE.nc
set -x

export XLFUNIT_31="sstoiqd.sst.$CDATE.grb"
export XLFUNIT_30="$MASK"
export XLFUNIT_34="$STDEV"
export XLFUNIT_51="${AMSR_name}day.grid.$CDATE.ieee"
export XLFUNIT_52="${AMSR_name}nte.grid.$CDATE.ieee"
export XLFUNIT_61="${AMSR_name}day.sobs.$CDATE.ieee"
export XLFUNIT_62="${AMSR_name}nte.sobs.$CDATE.ieee"


cat <<oiEOF > amsr.parm
sstoiqd.${AMSR_name}grid.$CDATE.nc
$iyrmid $immid $idmid
$cov_amsr
$cov_amsr2
oiEOF
cat amsr.parm

startmsg
$AMSRX <amsr.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating amsr superobs for $CDATE"
postmsg "$jlogfile" "$msg"

else
  msg="Warning - $AMSR_name file sstoiqd.${AMSR_name}grid.$CDATE.nc not available.  Skip."
  postmsg "$jlogfile" "$msg"
  continue
fi
 

done   # end daily loop generating insitu and satellite superobs.

if [ "$SENDCOM" = 'YES' ]; then
  [ -s buoy.sobs.$PDYm1.ieee ] && cp buoy.sobs.$PDYm1.ieee $COMOUT/sstoiqd.buoy.sobs.$PDYm1.ieee
  [ -s ship.cobs.$PDYm1.ieee ] && cp ship.cobs.$PDYm1.ieee $COMOUT/sstoiqd.ship.cobs.$PDYm1.ieee
fi


set +x
echo " "
echo "############# generate ICE superobs for $PDYm1. ########################"
set -x

kdays=0   # counter for number of days of ice data found
[ -s icecon.aggregate.ieee ] && rm icecon.aggregate.ieee

[ -s ice.dates ] && rm ice.dates
for CDATE in $ICE_DATElne_minus; do
  echo $CDATE >> ice.dates
done

for CDATE in `sort ice.dates`; do

iyrmid=`echo $CDATE | cut -c 1-4`
immid=`echo $CDATE | cut -c 5-6`
idmid=`echo $CDATE | cut -c 7-8`

file_in=sstoiqd.ice.$CDATE.grb

pgm=sstoiqd_ice2_1440x720
. prep_step

export XLFUNIT_11="sstoiqd.ice.$CDATE.grb"
export XLFUNIT_30="$MASK"
export XLFUNIT_51="ice.720x360.$CDATE.ieee"
export XLFUNIT_52="ice.1440x720.$CDATE.ieee"

cat << oiEOF > ice_read.parm
sstoiqd.ice.$CDATE.grb
$iyrmid $immid $idmid
oiEOF

startmsg
$ICEGRDX <ice_read.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating ice grid for $CDATE"
postmsg "$jlogfile" "$msg"

# add to aggregate file
if [ -s ice.1440x720.$CDATE.ieee ]; then
  cat ice.1440x720.$CDATE.ieee >> icecon.aggregate.ieee
  kdays=`expr $kdays + 1`
fi

done   # end daily loop generating ice grids


# using the grid fields generated above, compute median field and generate simulated ssts.

pgm=sstoiqd_ice2sst_med
. prep_step

export XLFUNIT_31="icecon.aggregate.ieee"
export XLFUNIT_20="$ICECOEF"
export XLFUNIT_30="$ICEMASK"
export XLFUNIT_40="$MASK"

export XLFUNIT_61="icesst.grid.ieee"
export XLFUNIT_71="ice.sobs.$PDYm1.ieee"
export XLFUNIT_81="icecon_med.grid.ieee"

cat << sstEOF > ice2sst.parm
icecon.aggregate.ieee
$iyrmid $immid $idmid $kdays
$cov_ice
$cov_ice2
sstEOF

startmsg
$ICE2SSTX <ice2sst.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating ice superobs"
postmsg "$jlogfile" "$msg"

if [ "$SENDCOM" = 'YES' ]; then
  [ -s ice.sobs.$PDYm1.ieee ] && cp ice.sobs.$PDYm1.ieee $COMOUT/sstoiqd.ice.sobs.$PDYm1.ieee
fi


set +x
echo " "
echo "############# SATELLITE BIAS CORRECTION. ########################"
set -x

#
# Initial data for EOT weights
minz=750       		# minimum number of smoothed zonal points
njj=10			# number of latitude points for zonal smoothing  
minob=1	        	# minimum number of obs used in zonal average
ioff=1			# number of longitude boxes smoothed for zonal sub
joff=1   		# number of latitude boxes smoothed for zonal sub
ndays=$ndays_satbc      # number of days used
drej=5.			# data REJECT LIMIT for ABSOLUTE SST increment 
crit=0.25	        # critical data/total area ratio 
modes=130	        # number of modes used <131
nby=7			# number of ship obs = 1 buoy obs 
nsh=1			# number of buoy obs = 1 
nst=1			# number of sat obs = 1 

#
# Initial data for EOT corrections
nuse=1                  # number of weights used (1, 3, or 5, default is 1) 
dfw1=1.0                # factor to smooth the weights: range 0 to 1    

#
# begin EOT weights for bias correction day & night
#
#  The 1st eotwt4 input record is 5 minz,njj,minob,ioff,joff
#  The 2nd eotwt4 input record is 5 parameters: nby,nsh,nst,ndays,nmodes
#  The 3rd record is 3 parameters: yrbuoy yrship yrsat
#  The 4th record is 3 parameters: daybuoy dayship daysat
#  The 5th record is 3 paramters: drej crit
#  The 6th input record is the date
#  All remaining records input filename templates which must be in order
#  The 5th file name is the input buoy data 
#  The 6th file name is the input ship data 
#  The 7th file name is the input satellite data
#  Multitple dates are read in and the complete file names
#    for these last 3 files are computed in the fortran code
#
#  Note all satellite fields, weights differ for day and night
#
#echo " "
#echo "minz njj minob ioff joff"
#echo "nby nsh nst ndays modes"
#echo "yrbuoy yrship yrsat"
#echo "daybuoy dayship daysat"
#echo "drej crit"
#echo "$PDYm1"
#echo " "


> sstoiqd.sat.list   # create empty file

BUOY=buoy.sobs.YYYYMMDD.ieee
SHIP=ship.cobs.YYYYMMDD.ieee

for SATname in $AVHRR_sats $AMSR_name; do
for typ in day nte; do

sattyp=$SATname$typ
SAT=$sattyp.sobs.YYYYMMDD.ieee
SATwts=$sattyp.wts.$PDYm1.ieee
SATtocorrect=`echo $SAT | sed "s/YYYYMMDD/${PDYm1}/g"`
if [ ! -s $SATtocorrect ];then
  msg="No file $SATtocorrect to correct.  Skip computation of correction for this $sattyp."
  postmsg "$jlogfile" "$msg"
  continue
fi

daybuoy=`echo $BUOY | awk '{print index($1,"YYYYMMDD")}'` # character count of buoy first number of date
dayship=`echo $SHIP | awk '{print index($1,"YYYYMMDD")}'` # character count of ship first number of date
daysat=`echo $SAT | awk '{print index($1,"YYYYMMDD")}'` # character count of sat first number of date
# not using yearly directories here. to minimize code changes for now, let eotwt4 apply year at date location as well.
yrbuoy=$daybuoy         # character count of buoy first number of year
yrship=$dayship         # character count of ship first number of year
yrsat=$daysat           # character count of sat first number of year

cat <<inxEOF> eotwt4-$sattyp.parm
$minz $njj $minob $ioff $joff
$nby $nsh $nst $ndays $modes
$yrbuoy $yrship $yrsat
$daybuoy $dayship $daysat
$drej $crit
$PDYm1
$BUOY
$SHIP
$SAT
inxEOF
cat eotwt4-$sattyp.parm

pgm=sstoiqd_eotbiaswt4
. prep_step
export XLFUNIT_32="$MASK2D"
export XLFUNIT_33="$CLIM2D"
export XLFUNIT_34="$FMODE"
export XLFUNIT_53="$SATwts"
export XLFUNIT_11="$BUOY"    # $BUOY is a template.  unit 11 gets opened for $ndays in executable 
export XLFUNIT_12="$SHIP"    # $SHIP is a template.  unit 12 gets opened for $ndays in executable 
export XLFUNIT_13="$SAT"     # $SAT is a template.  unit 13 gets opened for $ndays in executable 

startmsg
$SATBWTX <eotwt4-$sattyp.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done generating eot bias wts"
postmsg "$jlogfile" "$msg"


## apply the bias correction.

#
# begin EOT bias correction day & night
#
#
#  The 1st eotbias-cor4 input record is 2 parameter: nmodes, nuse
#  The 2nd eotbias-cor4 input record up nuse (up to 5 parameters): dfwN
#  The third record is the date to be run
#  The remaining record(s) is/are input wts.
#
#  Note all satellite fields, weights and biases differ for day and night

SATwts1=$sattyp.wts.$PDYm1.ieee
SATcobs=$sattyp.cobs.$PDYm1.ieee
ANALwts=$sattyp.awts.$PDYm1.ieee

cat <<inxEOF> eotcor4-$sattyp.parm
$modes $nuse
$dfw1
$PDYm1
$SATwts1
inxEOF
cat eotcor4-$sattyp.parm

pgm=sstoiqd_eotbiascor4
. prep_step

 echo  XLFUNIT_15="$SATwts1"  # unit 15 is opened in executable.  May be reused multiple files as determined by $nuse.
export XLFUNIT_13="$SATtocorrect"
export XLFUNIT_32="$MASK2D"
export XLFUNIT_34="$FMODE"
export XLFUNIT_39="$VARMD"

export XLFUNIT_51="$ANALwts"
export XLFUNIT_52="$SATcobs"

startmsg
$SATBCX <eotcor4-$sattyp.parm  >> $pgmout 2>errfile
export err=$?;err_chk
msg="Done applying bias correction to $sattyp"
postmsg "$jlogfile" "$msg"


echo sstoiqd.$ANALwts >> sstoiqd.sat.list  # passed in to analysis job
echo sstoiqd.$SATcobs >> sstoiqd.sat.list  # passed in to analysis job

if [ "$SENDCOM" = 'YES' ]; then
  if [ -s $ANALwts ]; then 
     cp $ANALwts $COMOUT/sstoiqd.$ANALwts
  else
     > $COMOUT/sstoiqd.$ANALwts  # create empty file
  fi
  if [ -s $SATcobs ];then
    cp $SATcobs $COMOUT/sstoiqd.$SATcobs
  else
    > $COMOUT/sstoiqd.$SATcobs   # create empty file
  fi
fi

done   # end typ (day/nte) loop
done   # end SATname loop

### end satellite bias correction

[ "$SENDCOM" = 'YES' ] && cp sstoiqd.sat.list $COMOUT

# restrict access to files containing ship ids
chgrp rstprod $COMOUT/sstoiqd.ships.*
chmod 640 $COMOUT/sstoiqd.ships.*

msg='ENDED NORMALLY.'
postmsg "$jlogfile" "$msg"

################## END OF SCRIPT #######################