#!/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=38
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 II: decode, calculate and write to grib 
# for the SPC SREF ensemble products (AWIPS)
#. prep_step

#startmsg
cp $PARMsref/sref_product_variable_SPC.tbl.$freq variable.tbl

# Calculate SPC's combined severe storm products
#get  flag/mask data 
#cp $SPCfix/spc_hicape_model_flag_gem hicape_model_flag_gem.out   #don't use land mask so as to get cptp in all areas

#get calibration data for cptp dryt rgn3

cp $SPCfix/spc_combine_pops.out combine_pops.out
cp $SPCfix/spc_combine_pops_dryt.out combine_pops_dryt.out
for x in 01 02 03 04 05 06 07 08 09 10 11 12 13 14 15 16 17 18 19 20 21 ; do
  cp $SPCfix/spc_combine_probs_svr_layer${x}.out combine_probs_svr_layer${x}.out
done

typeset -Z3 t0
typeset -Z3 t1
typeset -Z3 t2
typeset -Z2 h0
typeset -Z2 h1
typeset -Z2 h2

for rgn in _rgn3_ ; do
for hr in  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 ; do
#if [ $hr -le 36 ] ; then
 if [ $hr -le 39 ] ; then
      t0=$hr
      t1=`expr $hr - 1`
      t2=`expr $hr - 2`
      h0=$hr
      h1=`expr $hr - 1`
      h2=`expr $hr - 2`
      cp $SPCfix/spc_combine_pops${rgn}${cyc}_f${t0}.out.gz combine_pops${rgn}${cyc}_f${h0}.out.gz
      cp $SPCfix/spc_combine_pops${rgn}${cyc}_f${t0}.out.gz combine_pops${rgn}${cyc}_f${h1}.out.gz
      cp $SPCfix/spc_combine_pops${rgn}${cyc}_f${t0}.out.gz combine_pops${rgn}${cyc}_f${h2}.out.gz
      gunzip combine_pops${rgn}${cyc}_f${h0}.out.gz
      gunzip combine_pops${rgn}${cyc}_f${h1}.out.gz
      gunzip combine_pops${rgn}${cyc}_f${h2}.out.gz
 else
      t0=$hr
      h0=$hr
      cp $SPCfix/spc_combine_pops${rgn}${cyc}_f${t0}.out.gz combine_pops${rgn}${cyc}_f${h0}.out.gz
      gunzip combine_pops${rgn}${cyc}_f${h0}.out.gz
 fi
done
done
      
# Run product generator
cp filename.cluster0 filename
$EXECsref/sref_ens_gen < input >> $pgmout 2>errfile
export err=$?;err_chk

###################################################
# Part III: output grib data with proper file names
# (SPC one)
###################################################
if [ $freq = "3hrly" ]; then
 HR=00
else
 HR=01
fi
while [ $HR -le $ENDHOUR ]
do
 cp prob.sref.f$HR       $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob.f$HR
 cat prob.sref.f$HR >>   $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq

 rm -f mean.sref.f$HR
 rm -f spread.sref.f$HR
 rm -f prob.sref.f$HR
 rm -f max.sref.f$HR
 rm -f min.sref.f$HR
 rm -f mod.sref.f$HR
 rm -f p10.sref.f$HR
 rm -f p25.sref.f$HR
 rm -f p50.sref.f$HR
 rm -f p75.sref.f$HR
 rm -f p90.sref.f$HR

 HR=`expr $HR + $INCHOUR`
done

# Convert to grib2 format
CNVGRIB=${CNVGRIB:-/nwprod/util/exec/cnvgrib}
WGRIB2=${WGRIB2:-/nwprod/util/exec/wgrib2}
$CNVGRIB -g12 -p40 -nv $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2
$WGRIB2 $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2 -s >     $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2.idx

#if test "$SENDCOM" = "YES"
#then
#   mv sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq $COMOUT/.
#   if [ $SENDCOM_GB2 = YES ]
#   then
#   mv sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2 $COMOUT/.
#   mv sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2.idx $COMOUT/.
#   fi
#fi

if test "$SENDDBN" = "YES"
then
   $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB $job $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq

   #if [ $SENDDBN_GB2 = YES ]
   #then

   $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB_GB2 $job $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2

   $DBNROOT/bin/dbn_alert MODEL SREF_ENS_PGB_GB2_WIDX  $job $COMOUT/sref.t${cyc}z.pgrb${grid}_SPC.prob_$freq.grib2.idx

   #fi
fi

exit