#!/bin/ksh
#
#####################################################
# produce commom ensemble products (mean, spread and
# probability) of selected variables and write them 
# out in grib format
#
# programmer: Jun Du
# 07/10/2001: initial program
# 04/04/2002: J. Du, output 104 grid to ftp site
# 08/07/2003: J. Du, (1) 5 kfeta members added so products
#               is now 15 member based rather than 10 member 
#               based; (2) source code modified to follow 
#               NCEP ensemble grib extension standard
# 10/23/2003: J. Du, SREF member structure changed from
#               originally 3 5-member components (Eta,
#               kf_eta and RSM) to 2 component system
#               (10 Eta members and 5 RSM members)
# 02/22/2005: J. Du, added an extra program to produce similar
#                ensemble products -- mean, spread and prob --
#                but with more variables and levels as well as
#                outputs in individual forecast hours (instead of 
#                one big file) for the use by FVS (note: spread and
#                prob are produced but not saved since they are not 
#                used by FVS).
# 06/21/2005: J. Du, generalized the code to take any grids rather 
#                than one particular grid (212)
#
# 07/29/2005: B. Zhou, rewrite PART I 
#
# 10/26/2005: J. Du, add 6 WRF members and implement it in operational SREF
# 05/09/2006: J. Du, since different forecast region might have different 
#                requirement toward SREF products, separate Tables are used
#                for different grids (212, 216, 243)
# 06/29/2006: J. Du, add 2m-Td field to the fvs-version of SREF mean files
# 07/28/2008: J. Du, (1) increase WRF membership from 6 to 10 and reduce Eta
#                from 10 to 6; (2) add aviation ensemble products in standard
#                SREF products; and (3) add DTRA wind variance as a separate 
#                ensemble product
# 02/01/2011: J. Du and B. Zhou, added hourly output
# 04/01/2011: J. Du and B. Zhou, added 10-25-50-75-90%, mod and min/max
# 08/01/2011: B. Zhou, added SPC severe storm products
# 10/25/2011: J. Du and B. Zhou, added cluster mean
# 10/26/2011: J. Du, modified for the new SREF (NMMB, NMM, ARW)
# Hourly products: precip accumulation info is problematic for 2, 4, 6 ...hrs!
#####################################################
set -x

freq=$1
workdir=$3
method=NCEP
clusterdata=${COM_IN}/sref.$PDY/$cyc/cluster_$method

cd $workdir
######################################
# Part I: Fetching SREF data from /com 
######################################
# clean directory before job starts
rm -f r_gribawips.ctl*
rm -f r_gribawips.p01*
rm -f r_gribawips.p02*
rm -f r_gribawips.p03*
rm -f r_gribawips.p04*
rm -f r_gribawips.p05*
rm -f r_gribawips.p06*
rm -f r_gribawips.p07*
rm -f r_gribawips.p08*
rm -f r_gribawips.p09*
rm -f r_gribawips.p10*
rm -f r_gribawips.n01*
rm -f r_gribawips.n02*
rm -f r_gribawips.n03*
rm -f r_gribawips.n04*
rm -f r_gribawips.n05*
rm -f r_gribawips.n06*
rm -f r_gribawips.n07*
rm -f r_gribawips.n08*
rm -f r_gribawips.n09*
rm -f r_gribawips.n10*
rm -f *

export XLFRTEOPTS="namelist=old"

sh $utilscript/setup.sh

grid=$2
#grid=212
TOTMBR=21
KM=39             #vertical pressure levels in input grib files
yy=`echo $PDY | cut -c 1-4`
mm=`echo $PDY | cut -c 5-6`
dd=`echo $PDY | cut -c 7-8`

system=$4

if [ $freq = "3hrly" ]; then
ENDHOUR=87
INCHOUR=3
fi
if [ $freq = "1hrly" ]; then
ENDHOUR=39
INCHOUR=1
fi

if [ $freq = "3hrly" ]; then
fhrs="00 03 06 09 12 15 18 21 24 27 30 33 36 39 42 45 48 51 54 57 60 63 66 69 72 75 78 81 84 87"
interval=3
timestep=30
fi

if [ $freq = "1hrly" ]; then
fhrs="01 02 04 05 07 08 10 11 13 14 16 17 19 20 22 23 25 26 28 29 31 32 34 35 37 38"
interval=1
timestep=26
fi

# calculation based on the following order
# nmb_ctl: p01 
# nmb_n1:  n01
# nmb_p1:  p02
# nmb_n2:  n02
# nmb_p2:  p03
# nmb_n3:  n03  
# nmb_p3:  p04

# nmm_ctl: n04
# nmm_n1:  p05
# nmm_p1:  n05
# nmm_n2:  p06
# nmm_p2:  n06
# nmm_n3:  p07
# nmm_p3:  n07

# em_ctl:  p08
# em_n1:   n08
# em_p1:   p09
# em_n2:   n09
# em_p2:   p10
# em_n3:   n10
# em_p3:   ctl

MODEL[0]=
MODEL_PAIR[0]=
PAIR[0]=

MODEL[1]=sref_nmb
MODEL_PAIR[1]=ctl
PAIR[1]=p01
MODEL[2]=sref_nmb
MODEL_PAIR[2]=n1
PAIR[2]=n01

MODEL[3]=sref_nmb
MODEL_PAIR[3]=p1
PAIR[3]=p02
MODEL[4]=sref_nmb
MODEL_PAIR[4]=n2
PAIR[4]=n02

MODEL[5]=sref_nmb
MODEL_PAIR[5]=p2
PAIR[5]=p03
MODEL[6]=sref_nmb
MODEL_PAIR[6]=n3
PAIR[6]=n03

MODEL[7]=sref_nmb
MODEL_PAIR[7]=p3
PAIR[7]=p04
MODEL[8]=sref_nmm
MODEL_PAIR[8]=ctl
PAIR[8]=n04

MODEL[9]=sref_nmm
MODEL_PAIR[9]=n1
PAIR[9]=p05
MODEL[10]=sref_nmm
MODEL_PAIR[10]=p1
PAIR[10]=n05

MODEL[11]=sref_nmm
MODEL_PAIR[11]=n2
PAIR[11]=p06
#nmm members
MODEL[12]=sref_nmm
MODEL_PAIR[12]=p2
PAIR[12]=n06

MODEL[13]=sref_nmm
MODEL_PAIR[13]=n3
PAIR[13]=p07
MODEL[14]=sref_nmm
MODEL_PAIR[14]=p3
PAIR[14]=n07

MODEL[15]=sref_em
MODEL_PAIR[15]=ctl
PAIR[15]=p08
MODEL[16]=sref_em
MODEL_PAIR[16]=n1
PAIR[16]=n08

MODEL[17]=sref_em
MODEL_PAIR[17]=p1
PAIR[17]=p09
MODEL[18]=sref_em
MODEL_PAIR[18]=n2
PAIR[18]=n09

MODEL[19]=sref_em
MODEL_PAIR[19]=p2
PAIR[19]=p10
MODEL[20]=sref_em
MODEL_PAIR[20]=n3
PAIR[20]=n10

MODEL[21]=sref_em
MODEL_PAIR[21]=p3
PAIR[21]=ctl

MODEL[22]=
MODEL_PAIR[22]=
PAIR[22]=
#--------------------

typeset -Z2 HR
HR=00

NENS=0                               #total number of available members counted
nf=1
while [ $nf -le $TOTMBR ]            #loop for all members
do

#   cklast="no"
#   waittime=0
#   hasfile="yes"
#   while [ $cklast = "no" ] && [ $hasfile = "yes" ]                            # check if last fsct file is finished yet?
#   do
#    if [ -s $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$ENDHOUR ] && [ ${fsize[0]} -gt 800 ] ; then # last fsct file is finished?
#      cklast="yes"
#      #echo $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$ENDHOUR " exist"
#    else
#      if [ $waittime -gt 10 ] ; then                                         # over waiting time, last fcst file definitely not exist
#        hasfile="no"
#      fi                                                                      # last fcst file could come but waiting
#      echo "waiting " $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$ENDHOUR .....
#      sleep 60
#      waittime=`expr $waittime + 1`
#    fi
#   done

   HR=00
   hasfile="yes"
   while [ $HR -le $ENDHOUR ]
   do
    waittime=0
    cklast="no"
    while [ $cklast = "no" ] && [ $hasfile = "yes" ]  #if a fcst file at any single fcst time is missing, this member won't be counted
    do
     str=`ls -s $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR`  #get file size in block (1024 bytes)
     set -A fsize $str
     if [ -s ${COMIN}/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR ] && [ ${fsize[0]} -gt 800 ] ; then # fcst file is finished?
#      echo  ${COMIN}/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR " exist"
       cklast="yes"
     else
       sync
       if [ $waittime -gt 10 ] ; then                                       
         echo ${COMIN}/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR " is temporarily missing and will check one more time later!"
       fi                                                                   
       echo "waiting " ${COMIN}/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR .....
       sleep 60
       waittime=`expr $waittime + 1`
     fi
    done
    HR=`expr $HR + $INCHOUR`
   done
 
    HR=00
    while [ $HR -le $ENDHOUR ] && [ $hasfile = "yes" ]   #further check each fcst time files
    do
     ckall="no"
     icnt=1
     while [ $icnt -lt 1000 ]
     do
      str=`ls -s $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR`  #get file size in block (1024 bytes)
      set -A fsize $str
      if [ -s $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR ] && [ ${fsize[0]} -gt 800 ] ; then
       dump=0
       break
      else
#### add by WANG 20090212
       sync
       icnt=$((icnt + 1))
       sleep 5
      fi
      if [ $icnt -ge 120 ]
       then
       echo $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR "not exist"
       hasfile="no"
       msg="ABORTING after 10 minutes of waiting for $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR"
       err_exit $msg
      fi
     done
#### finish adding setion

     HR=`expr $HR + $INCHOUR`
    done
    echo ${MODEL[$nf]} ${MODEL_PAIR[$nf]} " done"  

    HR=00 
    if [ $hasfile = "yes" ] ; then
      while [ $HR -le $ENDHOUR ]
      do
        ln -s $COMIN/${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]}.f$HR r_gribawips.${PAIR[$nf]}.f$HR
        HR=`expr $HR + $INCHOUR`
      done
      member[$nf]=1
      NENS=`expr $NENS + 1`
    else
      member[$nf]=0
    fi

  nf=`expr $nf + 1`

done

echo $NENS $grid $KM $timestep $interval -1 > filename.cluster0
echo $fhrs >> filename.cluster0

nf=1
while [ $nf -le $TOTMBR ]
do
 if [ ${member[$nf]} -eq 1 ] ; then
   echo r_gribawips.${PAIR[$nf]} "  "  ${MODEL[$nf]}  "   -> " ${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]} >> filename.cluster0
 fi
 nf=`expr $nf + 1`
done

cat <<paramEOF >input
 &namin
   iyr=${yy},imon=${mm},idy=${dd},ihr=${cyc},ihh1=1,ihh2=${timestep},sys=$system
 &end
paramEOF

###############################################
# Part III: decode, calculate and write to grib
# for FVS verification system
###############################################
export pgm=sref_com_grib_fvs.$grid
#. prep_step

#startmsg
rm -f variable.tbl
cp filename.cluster0 filename
cp $PARMsref/sref_product_variable.tbl_fvs variable.tbl
$EXECsref/sref_ens_gen < input >> $pgmout 2>errfile
export err=$?;err_chk

HR=00
while [ $HR -le $ENDHOUR ]
do

#if test "$SENDCOM" = "YES"
#then
   mv mean.sref.f$HR   $COMOUT/sref.t${cyc}z.pgrb${grid}.mean_fvs.f$HR
#  mv spread.sref.f$HR $COMOUT/sref.t${cyc}z.pgrb${grid}.spread_fvs.f$HR
#  mv prob.sref.f$HR   $COMOUT/sref.t${cyc}z.pgrb${grid}.prob_fvs.f$HR
   rm spread.sref.f$HR
   rm prob.sref.f$HR
   rm max.sref.f$HR
   rm min.sref.f$HR
   rm mod.sref.f$HR
   rm p10.sref.f$HR
   rm p25.sref.f$HR
   rm p50.sref.f$HR
   rm p75.sref.f$HR
   rm p90.sref.f$HR
#fi

 if test "$SENDDBN" = "YES"
 then
   echo "This product will not be sent outside of the CCS"
#  $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB $job $COMOUT/sref.t${cyc}z.pgrb${grid}.mean_fvs.f$HR
#  $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB $job $COMOUT/sref.t${cyc}z.pgrb${grid}.spread_fvs.f$HR
#  $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB $job $COMOUT/sref.t${cyc}z.pgrb${grid}.prob_fvs.f$HR
 fi

HR=`expr $HR + $INCHOUR`
done

###############################################
# Part IV: wind variance product for DTRA project
###############################################
rm -f variable.tbl
cp filename.cluster0 filename
cp $PARMsref/sref_product_variable.tbl_DTRA variable.tbl
$EXECsref/sref_ens_gen < input >> $pgmout 2>errfile
export err=$?; err_chk

#if test "$SENDCOM" = "YES"; then
 HR=00
 while [ $HR -le $ENDHOUR ]
 do
  mv mean.sref.f$HR   $COMOUT/sref.t${cyc}z.pgrb${grid}.mean_DTRA.f$HR
  mv spread.sref.f$HR $COMOUT/sref.t${cyc}z.pgrb${grid}.spread_DTRA.f$HR
  rm prob.sref.f$HR
  rm max.sref.f$HR
  rm min.sref.f$HR
  rm mod.sref.f$HR
  rm p10.sref.f$HR
  rm p25.sref.f$HR
  rm p50.sref.f$HR
  rm p75.sref.f$HR
  rm p90.sref.f$HR
 HR=`expr $HR + $INCHOUR`
 done
#fi

# Release the BIASESTIMATE job
  $SMSBIN/setev release_qpfbiasestimate
exit

################################################################################
#          Ensemble Cluster Computation (submitted as a separate job!!!!)
################################################################################
if [ $freq = '3hrly' ]; then

rm -f cluster.*
rm -f variable.tbl
cp $PARMsref/sref_product_variable_cluster.tbl.212  variable.tbl
cp $clusterdata/sref_cluster_info sref_cluster_info

exec 2<&0 < sref_cluster_info

nfhr=1
for fhr in $fhrs ; do

 finished='no'
 while [ $finished = 'no' ] ; do

 read LINE
 set -A x $LINE

 if [ "${x[0]}" = $fhr ] ; then

 read LINE              #skip "**********"

  cluster=1
   while [ $cluster -le 6 ] ; do

    read LINE
    set -A y $LINE
    
    if [ ${y[4]} -gt 0 ] ; then
      echo ${y[4]} $grid $KM $timestep $interval $nfhr >> cluster.$cluster.f$fhr
      echo $fhrs >> cluster.$cluster.f$fhr
    
      n=1
      while [ $n -le ${y[4]} ] ; do
       read LINE  
       set -A z $LINE
       nf=${z[2]}
       echo r_gribawips.${PAIR[$nf]} "  "  ${MODEL[$nf]}  "   -> " ${MODEL[$nf]}.t${cyc}z.pgrb${grid}.${MODEL_PAIR[$nf]} >> cluster.$cluster.f$fhr  "mbr:" $nf
       n=`expr $n + 1 `
      done

      cp cluster.$cluster.f$fhr filename

      $EXECsref/sref_ens_gen < input >> output 2>errfile
      export err=$?; err_chk

      mv mean.sref.f$fhr $clusterdata/mean.sref.cluster$cluster.f$fhr
      mv spread.sref.f$fhr $clusterdata/spread.sref.cluster$cluster.f$fhr
      mv prob.sref.f$fhr $clusterdata/prob.sref.cluster$cluster.f$fhr
      rm -f max.sref.f$fhr min.sref.f$fhr mod.sref.f$fhr p10.sref.f$fhr p25.sref.f$fhr p50.sref.f$fhr p75.sref.f$fhr p90.sref.f$fhr

      echo "Forecast hour " $fhr "-> Cluster " $cluster         
    else
      echo "cluster" $cluster " has no forecasts"
    fi 

      cluster=`expr $cluster + 1 `

  done

  finished='yes'

 fi 

 done

echo    "Forecast hout " $fhr "done ..."

nfhr=`expr $nfhr + 1 `

done

fi

exit