#!/bin/ksh
#
# LSF batch script to run an MPI application
#
#BSUB -P P48503002             # project code
#BSUB -W 00:10                 # wall-clock time (hrs:mins)
#BSUB -n 16                     # number of tasks in job         
#BSUB -R "span[ptile=16]"      # run 16 MPI tasks per node
#BSUB -J gsi                   # job name
#BSUB -o gsi.%J.out          # output file name in which %J is replaced by the job ID
#BSUB -e gsi.%J.err          # error file name in which %J is replaced by the job ID
## #BSUB -q regular               # queue
#BSUB -q small               # queue

#

set -x

ROOTDIR=/glade/p/work/mhu/gsi/rtma/rtma/RTMA
work_dir=/glade/p/work/mhu/gsi/rtma/rtma/RTMA/rpostprd
fixparm=${ROOTDIR}/fix
rtmagsidir=$ROOTDIR/rtmaprd
rtmafgdir=$ROOTDIR/rtmagus

CYCLE=2012052811
RUN_COMMAND="mpirun.lsf"

#========================================================================
lstats=$work_dir/lstats
tmpdir=${rtmagsidir}
. ${ROOTDIR}/rtma_envir.sh
uncompress=gunzip

  YYYYMMDD=`echo ${CYCLE} | cut -c 1-8`
  HH=`echo ${CYCLE}   | cut -c 9-10`
  cycle_date="${YYYYMMDD} ${HH}"

#========================================================================
# (1) read in from original unformatted gsi ob stats files and write out
# formatted, streamlined versions of these files. Besides (lat,lon),
# output files will also contain location in grid units (xx,yy) for each
# observation.
#
# (2) read in bckg std error files for psi,chi,ps,tv,pseudorh
# in (ny,nx) configuration and (1) transpose files to
# (nx,ny) configuration and (2) get estimates of error files
# for u,v,ps,t,q.
#
# (3) perform cressman or barnes analysis to get gridded estimate of
# analysis errors for u,v,t,ps,q
#========================================================================

 if [ ! -d $work_dir ] ; then
    mkdir -p $work_dir
 fi

 cd $work_dir

 if [ -e diag_conv_anl.$CYCLE ] ; then
      rm diag_conv_anl.$CYCLE
 fi
 cp  $tmpdir/diag_conv_anl.$CYCLE .

 cp  $tmpdir/bckg_dxdy.dat .
 cp  $tmpdir/bckg_qsat.dat .
 cp  $tmpdir/bckg_psfc.dat .
 cp  $tmpdir/bckgvar.dat_psi .
 cp  $tmpdir/bckgvar.dat_chi .
 cp  $tmpdir/bckgvar.dat_ps bckgvar.dat_ps0
 cp  $tmpdir/bckgvar.dat_t bckgvar.dat_t0
 cp  $tmpdir/bckgvar.dat_pseudorh .
 cp  $tmpdir/bckgvar.dat_gust bckgvar.dat_gust0
 cp  $tmpdir/bckgvar.dat_vis  bckgvar.dat_vis0
#cp  $tmpdir/bckgvar.dat_pblh bckgvar.dat_pblh0
#cp $tmpdir/bckgvar.dat_dist  bckgvar.dat_dist0

 cp  $tmpdir/bckg_z.dat .

 cp ${fixgsi}/errfield.dat_clim_2008042119 errfield.dat_clim
 ls -lrt $tmpdir/gradx.dat_01_* | grep -c gradx.dat > used_iterations.dat

 nfiles=`cat used_iterations.dat`
 ic=0
 while [ $ic -le $(($nfiles-1)) ] ; do
     fic=`printf %03d $ic`
     ln -s -f  $tmpdir/gradx.dat_01_${fic} gradx.dat_01_${fic}
     ln -s -f  $tmpdir/grady.dat_01_${fic} grady.dat_01_${fic}
     let "ic=ic+1"
 done
 ln -s -f  grady.dat_01_000  grady.dat

 cp $ROOTDIR/fix/rtma_${cgrid}_latlon_mpfactor_slmask.data rtma_latlon_mpfactor.dat

if [ ${anl_or_01} = 01 ] ; then
  FHH=01
fi

if [ ${anl_or_01} = 02 ] ; then
  FHH=02
fi

if [[ "${anl_or_01}" != "01" && "${anl_or_01}" != "02" ]] ; then
  FHH=00
fi

cat << EOF > ${work_dir}/bckgerr_infl_input
 &bckgerr_infl
   rinflu=1.15,
   rinflv=1.15,
   rinflt=1.0,
   rinflq=0.4,
   rinflp=1.0,
   itdiff=$FHH
   efold=1.5
/
EOF
#************************************************************
#
#************************************************************
#==> generate input file used for cressman/barnes analysis

cat << EOF > ${work_dir}/cress_barnesparm.anl
 &cress_barnes
   radius_u=20.,rvicinity_u=0.5,rinovmax_u=100.,
   radius_v=20.,rvicinity_v=0.5,rinovmax_v=100.,
   radius_t=18.,rvicinity_t=0.5,rinovmax_t=100.,
   radius_p=18.,rvicinity_p=0.5,rinovmax_p=10000.,
   radius_q=18.,rvicinity_q=0.5,rinovmax_q=100.
 /
 &cress_barnes_anlqlty
   cgrid='${cgrid}',npass=1,
   radius_u=12.,rinovmax_u=35.,
   radius_v=12.,rinovmax_v=35.,
   radius_t=12.,rinovmax_t=15.,
   radius_p=10.,rinovmax_p=3000.,
   radius_q=10.,rinovmax_q=500.,
   radius_gust=15.,rinovmax_gust=35.,
   radius_vis=20.,rinovmax_vis=200000.,
   lpadjust=.true.,
   usebckg=.true.,
   cvmodel(1)='sigges',
   cvmodel(2)='siganl',
   cvmodel(3)='sigcress',
   cvmodel(4)='siganlerr',
   rmusecv(1)=-3.,
 /
EOF

#==> generate input file used for lanczos analysis error

cat << EOF > ${work_dir}/lanczosparm.anl
 &anlerrmethod
   projway=.true.
   climatology=.false.,
   anlqlty=.false.,
/
EOF

ln -sf sigges siganlerr

#************************************************************
#
#************************************************************
#==> generate input file needed for conversion to grib format

if [ ${anl_or_01} = 01 ] ; then
  CYCLEm1=`date +%Y%m%d%H -d "${cycle_date} 1 hours ago"`
  YYYY=`echo ${CYCLEm1} | cut -c 1-4`
  MM=`echo ${CYCLEm1}   | cut -c 5-6`
  DD=`echo ${CYCLEm1}   | cut -c 7-8`
  HH=`echo ${CYCLEm1}   | cut -c 9-10`
  FHH=01
  LSTATUS=.true.
fi 

if [ ${anl_or_01} = 02 ] ; then
  CYCLEm2=`date +%Y%m%d%H -d "${cycle_date} 2 hours ago"`
  YYYY=`echo ${CYCLEm2} | cut -c 1-4`
  MM=`echo ${CYCLEm2}   | cut -c 5-6`
  DD=`echo ${CYCLEm2}   | cut -c 7-8`
  HH=`echo ${CYCLEm2}   | cut -c 9-10`
  FHH=02
  LSTATUS=.true.
fi

if [[ "${anl_or_01}" != "01" && "${anl_or_01}" != "02" ]] ; then
  YYYY=`echo ${CYCLE} | cut -c 1-4`
  MM=`echo ${CYCLE}   | cut -c 5-6`
  DD=`echo ${CYCLE}   | cut -c 7-8`
  HH=`echo ${CYCLE}   | cut -c 9-10`
  FHH=00
  LSTATUS=.false.
fi

usefgat=.false.
 
cat << EOF > parm_ptime_namelist
&ptimeinfo
    iyear01=${YYYY},
    imonth01=${MM},
    iday01=${DD},
    ihour01=${HH},
    ifcsthr=${FHH},
    guessfcst=$LSTATUS
    usefgat=$usefgat
/
EOF

cat << EOF > hgridoptions_input
&hgridoptions
    lhalfres=.true.,
    iopt=1,
    npass=1,
    iopthgt=0,
    npasshgt=0,
    smcf1=0.50,
    smcf2=0.52,
    one2one=.false.,
/
EOF
#************************************************************
#
#************************************************************
/bin/rm -rf fort.9 fort.19 fort.29
ln -s -f diag_conv_anl.$CYCLE diag_conv.dat
ln -s -f  parm_ptime_namelist fort.9
ln -s -f  parm_ptimem1_namelist fort.19
ln -s -f  parm_ptimep1_namelist fort.29
ln -s -f  $fixparm/ruc2_grib.table_V2 GRIBTABLE

if [ ! -e diag_conv_ges.$CYCLE ] ; then
  cp $tmpdir/diag_conv_ges.$CYCLE .
fi

if [ ! -e diag_conv_02.$CYCLE ] ; then
  cp  $tmpdir/diag_conv_02.$CYCLE .
fi

if [ ! -e diag_conv_03.$CYCLE ] ; then
  cp  $tmpdir/diag_conv_03.$CYCLE .
fi


ln -s -f diag_conv_ges.$CYCLE diag_conv_ges.dat
ln -s -f diag_conv_02.$CYCLE diag_conv_02.dat
ln -s -f diag_conv_anl.$CYCLE diag_conv_anl.dat


#Note:
# bmax(1) - maximum allowed bias for wind speed (m/s)
# bmax(2) - maximum allowed bias for temperature (K)
# bmax(3) - maximum allowed bias for psfc (Pa)
# bmax(4) - maximum allowed bias for spec. hum. (g/kg)
# bmax(5) - maximum allowed bias for wind speed (satellite data) (m/s)
# nt      - number of additional past cycles
#

  cyc1=`date +%Y%m%d%H -d "${cycle_date} 1 hours ago"`
  cyc2=`date +%Y%m%d%H -d "${cycle_date} 2 hours ago"`
  cyc3=`date +%Y%m%d%H -d "${cycle_date} 3 hours ago"`
  cyc4=`date +%Y%m%d%H -d "${cycle_date} 4 hours ago"`
  cyc5=`date +%Y%m%d%H -d "${cycle_date} 5 hours ago"`
  cyc6=`date +%Y%m%d%H -d "${cycle_date} 6 hours ago"`

  nt0=6

cat << EOF > ${work_dir}/rj_param_input
 &rj_param
   tcyc='${CYCLE}',
   wmax280=5.0,  tmax180=10., pmax180=3000., qmax180=5.e-03, spdmax280=15.,
   wmax281=5.0,  tmax181=10., pmax181=3000., qmax181=5.e-03, spdmax281=15.,
   wmax282=5.0,  tmax182=10., pmax182=3000., qmax182=5.e-03, spdmax282=15.,
   wmax283=5.0,  tmax183=10., pmax183=3000., qmax183=5.e-03, spdmax283=15.,
   wmax284=5.0,  tmax184=10., pmax184=3000., qmax184=5.e-03, spdmax284=15.,
   wmax285=5.0,  tmax185=10., pmax185=3000., qmax185=5.e-03, spdmax285=15.,
   wmax286=5.0,  tmax186=10., pmax186=3000., qmax186=5.e-03, spdmax286=15.,
   wmax287=100., tmax187=10., pmax187=3000., qmax187=5.e-03, spdmax287=15.,
   wmax288=3.0,  tmax188=10., pmax188=3000., qmax188=5.e-03, spdmax288=15.,
   wmax289=5.0,  tmax189=10., pmax189=3000., qmax189=5.e-03, spdmax289=15.,
   wmax290=5.0,  tmax190=10., pmax190=3000., qmax190=5.e-03, spdmax290=15.,
   wmax291=5.0,  tmax191=10., pmax191=3000., qmax191=5.e-03, spdmax291=15.,
   wmax292=5.0,  tmax192=10., pmax192=3000., qmax192=5.e-03, spdmax292=15.,
   wmax293=5.0,  tmax193=10., pmax193=3000., qmax193=5.e-03, spdmax293=15.,
   wmax294=5.0,  tmax194=10., pmax194=3000., qmax194=5.e-03, spdmax294=15.,
   wmax295=5.0,  tmax195=10., pmax195=3000., qmax195=5.e-03, spdmax294=15.,
   ntw=1,
   ntt=3,
   ntp=3,
   ntq=3,
   ntspd=3,
   cycall(1)='${cyc1}',
   cycall(2)='${cyc2}',
   cycall(3)='${cyc3}',
   cycall(4)='${cyc4}',
   cycall(5)='${cyc5}',
   cycall(6)='${cyc6}'
/
EOF

   if [ ! -d ${lstats} ] ; then
      mkdir -p ${lstats}
   fi
   ges_or_anl=anl

   if [ ! -e diag_conv_ges.$CYCLE ] ; then
     cp $tmpdir/diag_conv_ges.$CYCLE.Z .
     ${uncompress} diag_conv_ges.$CYCLE.Z
   fi

   if [ ! -e diag_conv_anl.$CYCLE ] ; then
     cp $tmpdir/diag_conv_anl.$CYCLE.Z .
     ${uncompress} diag_conv_anl.$CYCLE.Z
   fi

   ges_or_anl=ges

   cp diag_conv_${ges_or_anl}.${CYCLE} diag_conv_for_rj.dat

   
   it=1
   while [ $it -le $nt0 ] ; do
     if [ $it -eq 1 ] ; then cycit=$cyc1 ; fi ; 
     if [ $it -eq 2 ] ; then cycit=$cyc2 ; fi ; 
     if [ $it -eq 3 ] ; then cycit=$cyc3 ; fi ; 
     if [ $it -eq 4 ] ; then cycit=$cyc4 ; fi ; 
     if [ $it -eq 5 ] ; then cycit=$cyc5 ; fi ; 
     if [ $it -eq 6 ] ; then cycit=$cyc6 ; fi ; 
  
     if [ -e ${lstats}/stats_w_rj.out_${cycit} ] ; then
       cp ${lstats}/stats_w_rj.out_${cycit} .
     fi

     if [ -e ${lstats}/stats_t_rj.out_${cycit} ] ; then
       cp ${lstats}/stats_t_rj.out_${cycit} .
     fi

     if [ -e ${lstats}/stats_ps_rj.out_${cycit} ] ; then
       cp ${lstats}/stats_ps_rj.out_${cycit} .
     fi

     if [ -e ${lstats}/stats_q_rj.out_${cycit} ] ; then
       cp ${lstats}/stats_q_rj.out_${cycit} .
     fi
     
     if [ -e ${lstats}/stats_spd_rj.out_${cycit} ] ; then
       cp ${lstats}/stats_spd_rj.out_${cycit} .
     fi

     let "it=it+1"
   done


if [ -e $tmpdir/t_rejectlist ] ; then
  cp $tmpdir/t_rejectlist .
fi

if [ -e $tmpdir/q_rejectlist ] ; then
  cp $tmpdir/q_rejectlist .
fi

if [ -e $tmpdir/p_rejectlist ] ; then
  cp $tmpdir/p_rejectlist .
fi

if [ -e $tmpdir/w_rejectlist ] ; then
  cp $tmpdir/w_rejectlist .
fi

if [ -e $tmpdir/t_day_rejectlist ] ; then
  cp $tmpdir/t_day_rejectlist .
fi

if [ -e $tmpdir/t_night_rejectlist ] ; then
  cp $tmpdir/t_night_rejectlist .
fi

if [ -e $tmpdir/q_day_rejectlist ] ; then
  cp $tmpdir/q_day_rejectlist .
fi

if [ -e $tmpdir/q_night_rejectlist ] ; then
  cp $tmpdir/q_night_rejectlist .
fi

if [ -e post_done ] ; then
    rm post_done
fi

cp $ROOTDIR/rtma_post/rtma_post.exe .
cp $rtmafgdir/gridname_input .
cp $rtmafgdir/twodvar_input_bi sigges
cp $rtmafgdir/slabs2.dat .
cp $rtmafgdir/slabs2_nobiasc.dat .
cp $rtmagsidir/wrf_inout siganl

# run
${RUN_COMMAND} ${work_dir}/rtma_post.exe <  ${work_dir}/cress_barnesparm.anl > stdout_post

list="ps t q u v spd gust vis"
for type in $list; do
    mv ${type}_obs.listing_iter_01      ${type}_obs.listing_iter_01.${CYCLE}
    mv ${type}_obs.listing_iter_02      ${type}_obs.listing_iter_02.${CYCLE}
    mv ${type}_obs.listing_iter_anl     ${type}_obs.listing_iter_anl.${CYCLE}
    mv ${type}_obs.listing_iter_01_aux  ${type}_obs.listing_iter_01_aux.${CYCLE}
done


list="w t p q spd mass wind"
for type in $list; do
    cp ${type}_rjlist.txt_dynamic_$CYCLE     ${lstats}/.
done

list="w t ps q spd"
for type in $list; do
    cp stats_${type}_rj.out          ${lstats}/stats_${type}_rj.out_$CYCLE
done

if [[ "$lbiascor" = ".true." ]] ; then
   cp rtma_biascor_out.dat_$CYCLE    ${lstats}/rtma_biascor_out.dat_$CYCLE
fi

exit 0
