#!/usr/bin/csh
set echo on
date
setenv yyyymmdd `date +%Y%m%d`
setenv mtime `date +%H`
#if( $mtime >= 12 & $mtime <= 23 ) then
#setenv ens_hr  09
#else
#setenv ens_hr 21
#endif
#setenv yyyymmdd 20121124
setenv ens_hr  09
setenv JOBHOME "/meso/save/wx20jd/research/sensitivity_sref/work/"             ###Minghua Comment: modify to your own directory.
setenv JOBHOME_sed "\/meso\/save\/wx20jd\/research\/sensitivity_sref\/work\/"   ###Minghua Comment: modify to your own directory.
setenv SREFDIR  /com/sref/prod/sref.${yyyymmdd}/${ens_hr}/pgrb         ###Minghua Comment: modify to your own CMC data set directory.
#setenv SREFDIR  /com/sref/prod/sref.${yyyymmdd}/global_pgrb/${ens_hr}        ###Minghua Comment: modify to your own CMC data set directory.
setenv GET_BIN   $JOBHOME"GET_BIN"
setenv GET_BIN_sed "${JOBHOME_sed}GET_BIN\/"
setenv GET_EOF_sed "${JOBHOME_sed}EOF\/"
setenv VARIANCE $JOBHOME"VARIANCE"
setenv EOF $JOBHOME"EOF"
setenv SENSITIVITY $JOBHOME"SENSITIVITY"
setenv SENSITIVITY_sed "${JOBHOME_sed}SENSITIVITY\/"
setenv MAIN_CSH $JOBHOME"MAIN_CSH"
setenv OUTPUT $JOBHOME"OUTPUT"
setenv center 
setenv region                                                   ###Minghua Comment: region: "1"==>30N-50N, 95W-55W; "2"==>
setenv day
setenv ndate ${yyyymmdd}${ens_hr}
setenv PSL_SREF "psl_all_SREF.dat"
setenv Z500_SREF "z500_all_SREF.dat"
setenv endate  "/nwprod/util/exec/ndate"  ###Minghua Comment: this software is for computing valid time.
setenv srflast ${SREFDIR}/sref_nmb.t${ens_hr}z.pgrb132.p3.f84.grib2       ###Minghua Comment:check whether data available

###############################################################################
###################              Set Parameters            ####################
###############################################################################
setenv tt  "8"                        ###The sample data has 8 forecast times.
setenv tinc "12"                      ###The interval between two time steps
setenv icent "4"                      ###Minghua Comment: "1"==>NCEP,CMC and NCEP+CMC; "2"==>NCEP,CMC,NCEP+CMC and ECMWF; "3"==>NCEP; "4"==>SREF; "5"==>ECMWF; "6"==>CMC,NCEP+CMC  
setenv iday  "1"                      ###Minghua Comment: "1"==>6,5,4,3,2,1; "2"==>6,5,4; "3"==>3,2,1 
#setenv vlat1 "10"                    ###Minghua Comment: default starting lat for variance plots
setenv vlat1 "8"                      ###For SREF data 
#setenv vlat2 "70"
setenv vlat2 "78"                     ###Minghua Comment: default ending lat for variance plots
#setenv vlon1 "210"                   ###Minghua Comment: default staring lon for variance plots
setenv vlon1 "-167"                   ###For SREF data, must be negative.
#setenv vlon2 "320"                   ####Minghua Comment: default ending lon for variance plots   
setenv vlon2 "-47"                    ###Minghua Comment: Must be negative for SREF data.
setenv depdata "dep_MSLP"
setenv e1 "1"
setenv e2 "21"   ##for test I chose 7, the csh file in ./GET_BIN/ should be checked.
# reso=16km for SREF (how to set reso?)  ###Minghua Comment: GrADS automatically generates gridded data instead of using km.   
setenv reso "1"
#setenv reso "0.5"                 ###This means half degree lat/lon. I tested several reso and think this reso is good in terms of efficiency and information. 
###############################################################################
###################           Set center combination       ####################       
###############################################################################
if ($icent == 1) then
  setenv icstr "NCEP CMC NCEP+CMC"
  setenv datastr "-e $ncplast & -e $cmclast"
  setenv memn "20"
  setenv memc "20"
  setenv meme "0"
  setenv memall "40"
endif
if ($icent == 2) then
    setenv icstr "NCEP CMC NCEP+CMC ECMWF"
    setenv datastr "-e $ncplast & -e $cmclast & -e $eclast"
    setenv memn "20"
    setenv memc "20"
    setenv meme "51"
    setenv memall "40"
endif
if($icent == 3) then
      setenv icstr "NCEP"
      setenv datastr "-e $ncplast"
      setenv memn "20"
      setenv memc "0"
      setenv meme "0"
      setenv memall "20"
endif
if($icent == 4) then 
          setenv icstr "SREF"
          setenv datastr "-e $srflast"
          setenv memn "21"
          setenv memc "0"
          setenv meme "0"
          setenv memall "21"

endif
if($icent == 5) then
          setenv icstr "ECMWF"
          setenv datastr "-e $eclast"
          setenv memn "0"
          setenv memc "0"
          setenv meme "51"
          setenv memall "51"
endif
if($icent == 6) then
          setenv icstr "CMC NCEP+CMC"
          setenv datastr "-e $ncplast & -e $cmclast"
          setenv memn "20"
          setenv memc "20"
          setenv meme "0"
          setenv memall "40"
endif

if($iday == 1) then
#         setenv idstr "6 5 4 3 2 1"
          setenv idstr "3 2 1"
endif
###############################################################################
###################  Check whether data is there  #############################
###############################################################################
set nn=1
while ($nn <= 10)
  if($datastr) then

###############################################################################
################### Get binary data from single grib2 data  ###################
###############################################################################
cd $GET_BIN
rm *.dat
csh getbin_2times_SREF.csh

###############################################################################
#########    Compute variance at each time step and plot figures     ##########
###############################################################################
cd $GET_BIN
foreach center ($icstr) #bigin loop center
rm psl_all_${center}.ctl
rm z500_all_${center}.ctl
cat ${GET_BIN}/psl_all_${center}_samp.ctl | sed -e "s/psl_all_${center}.dat/"${GET_BIN_sed}"psl_all_${center}.dat/g" >file.temp && mv file.temp ${GET_BIN}/psl_all_${center}.ctl
cat ${GET_BIN}/z500_all_${center}_samp.ctl | sed -e "s/z500_all_${center}.dat/"${GET_BIN_sed}"z500_all_${center}.dat/g" >file.temp && mv file.temp ${GET_BIN}/z500_all_${center}.ctl
end  #end loop center

cd $GET_BIN
if ($icent == 1 | $icent == 2 |$icent == 6) then
/usrx/local/grads/bin/2.0.a3/xgrads -blc 'run combine.gs'
/usrx/local/grads/bin/2.0.a3/xgrads -blc 'run gen_ano.gs'
endif

foreach center ($icstr)                                         ### Begin loop 1 for center


if ($center == "NCEP" | $center == "SREF") then
  setenv e2 21                    ##for test I chose 7.
  else
  if ($center == "ECMWF") then
    setenv e2 51
  endif
  if ($center == "NCEP+CMC") then
    setenv e2 40
  endif
endif

cd $OUTPUT
mkdir $ndate
cd $ndate
cp ${JOBHOME}cbarn.gs  $OUTPUT/${ndate}/cbarn.gs
cp ${JOBHOME}color-var.gs  $OUTPUT/${ndate}/color-var.gs
/usrx/local/grads/bin/2.0.a3/xgrads -bpc "run $VARIANCE/var.gs $tt $vlat1 $vlat2 $vlon1 $vlon2 $memn $memc $meme $memall $GET_BIN/psl_all_${center}.ctl $tinc $e1 $e2 $reso $ndate $center" 
foreach region (1 2)                                            ### Begin loop2 for region 


if ($region == 1) then
setenv lat1 30
setenv lat2 50
setenv lon1 265
setenv lon2 295
    else
if ($region == 2) then
setenv lat1 32
setenv lat2 45
setenv lon1 281
setenv lon2 298
    endif
endif

#setenv xgrid `expr \( $lon2 - $lon1 \) / $reso + 1`  
setenv xgrid `echo $lon1 $lon2 $reso | awk ' { printf( "%3.0f\n",  ($2-$1)/$3+1)} '`
setenv ygrid `echo $lat1 $lat2 $reso | awk ' { printf( "%3.0f\n",  ($2-$1)/$3+1)} '`
setenv agrid `expr $xgrid \* $ygrid`
mkdir ${OUTPUT}/${ndate}/"region"$region

foreach day ($idstr)                                      ### Begin loop3 for valid day

mkdir  ${OUTPUT}/${ndate}/"region"$region/"day"${day}
setenv ttt `expr $day \* 2 + 1`
setenv fcsthr ` expr \( $ttt - 1 \) \* 12 `
setenv valid ` $endate $fcsthr $yyyymmdd$ens_hr`
cd ${OUTPUT}/${ndate}
cp ${JOBHOME}cbarn.gs  $OUTPUT/${ndate}/cbarn.gs
cp ${JOBHOME}color-var.gs  $OUTPUT/${ndate}/color-var.gs
/usrx/local/grads/bin/2.0.a3/xgrads -bpc "run $VARIANCE/dep_MSLP.gs $ttt $lat1 $lat2 $lon1 $lon2 $memn $memc $meme $memall $GET_BIN/psl_all_${center}.ctl $depdata $e1 $e2 $reso $ndate $center $valid "

mv Valid_*.gif ${OUTPUT}/${ndate}/"region"${region}/"day"${day}/


 
#######################################################################
###################Compute EOF Principal Components.###################
#######################################################################
mv $depdata".dat"  $EOF
cd $EOF
sed s/113/$e2/g param_sample.txt > param_temp.txt
sed s/6771/$agrid/g param_temp.txt > param.txt
rm param_temp.txt
cp EOF_output.for EOF_output.f
f77 EOF_output.f  ###Please modify this to your own fortran compiler
./a.out
rm a.out  EOF_output.f
rm $depdata".dat"

foreach inum (1 2 3)                                 
rm time${inum}norm.txt
rm space${inum}.txt
rm time${inum}.txt
rm yprime${inum}norm.txt
end


###############################################################################
##################    plot EOF patterns      ##################################
###############################################################################
foreach inum (1 2 3)
rm eof${inum}_${center}.ctl
sed s/__VARI/yprime${inum}norm.dat/g eof_template.ctl > eof${inum}_${center}.ctl
echo "XDEF $xgrid linear $lon1 $reso " >> eof${inum}_${center}.ctl
echo "YDEF $ygrid linear $lat1 $reso " >> eof${inum}_${center}.ctl
cat var.ctl >> eof${inum}_${center}.ctl
end

setenv eofvar1 "`sed -n 1p eof_var.txt`"
setenv eofvar2 "`sed -n 2p eof_var.txt`"
setenv eofvar3 "`sed -n 3p eof_var.txt`"
/usrx/local/grads/bin/2.0.a3/xgrads -bpc "run eof.gs ${center} $OUTPUT ${ndate} $xgrid $ygrid $valid $lon1 $lon2 $lat1 $lat2 $eofvar1 $eofvar2 $eofvar3" 
rm eof_var.txt
mv ${OUTPUT}/${ndate}/Valid_EOF*.gif ${OUTPUT}/${ndate}/"region"${region}/"day"${day}/

###############################################################################
#####################       Compute Sensitivity      ##########################
###############################################################################
cd $OUTPUT/
cd $ndate
foreach inum (1 2 3)                                            ### Begin loop 4 for three PCs.  
cp ${JOBHOME}cbarn.gs  $OUTPUT/${ndate}/cbarn.gs
if ( $e2 == "21" ) then                                         ###If for colorbar, only 95% significant value will be plot 
cp ${JOBHOME}color.gs  $OUTPUT/${ndate}/color.gs
   else
if ($e2 == "51" ) then
cp ${JOBHOME}color-ec.gs  $OUTPUT/${ndate}/color.gs
endif
if ($e2 == "40") then
cp  ${JOBHOME}color-all.gs  $OUTPUT/${ndate}/color.gs
endif  
endif                                                           ###Endif for colorbar

if ( $e2 == "21" | $e2 == "51"  ) then
  sed s/__VARI/"${GET_EOF_sed}"time${inum}norm.dat/g $JOBHOME"time_template.ctl" > $SENSITIVITY/time${inum}norm.ctl
  echo "XDEF 1 linear 275.0 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
  echo "YDEF 1 linear 30 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
  cat ${JOBHOME}/member.ctl >> $SENSITIVITY/time${inum}norm.ctl
  cp $SENSITIVITY/sen_template.gs  $SENSITIVITY/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/113/${e2}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/999/${ttt}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/timenorm.ctl/"${SENSITIVITY_sed}"time${inum}norm.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/psl_all.ctl/"${GET_BIN_sed}"psl_all_${center}.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${GET_BIN}/psl_all_${center}_samp.ctl | sed -e "s/psl_all_${center}.dat/"${GET_BIN_sed}"psl_all_${center}.dat/g" >file.temp && mv file.temp ${GET_BIN}/psl_all_${center}.ctl
  /usrx/local/grads/bin/2.0.a3/xgrads -blc "run ${SENSITIVITY}/sen_eof${inum}.gs $center $ndate $inum $valid $lon1 $lon2 $lat1 $lat2 "
  rm cbarn.gs
  rm color.gs
  rm ${SENSITIVITY}/time${inum}norm.ctl
  rm ${SENSITIVITY}/sen_eof${inum}.gs
  rm *.dat
else
  sed s/__VARI/"${GET_EOF_sed}"time${inum}norm.dat/g $JOBHOME"time_template.ctl" > $SENSITIVITY/time${inum}norm.ctl
  echo "XDEF 1 linear 275.0 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
  echo "YDEF 1 linear 30 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
  cat ${JOBHOME}/member.ctl >> $SENSITIVITY/time${inum}norm.ctl
  cp $SENSITIVITY/sen_template_all.gs  $SENSITIVITY/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/113/${e2}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/999/${ttt}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/timenorm.ctl/"${SENSITIVITY_sed}"time${inum}norm.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/psl_all.ctl/"${GET_BIN_sed}"psl_all_${center}.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${SENSITIVITY}/sen_eof${inum}.gs | sed -e "s/psl_ano.ctl/"${GET_BIN_sed}"psl_all_${center}_ano.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}.gs
  cat ${GET_BIN}/psl_all_${center}_ano_samp.ctl | sed -e "s/psl_all_${center}_ano.dat/"${GET_BIN_sed}"psl_all_${center}_ano.dat/g" >file.temp && mv file.temp ${GET_BIN}/psl_all_${center}_ano.ctl
  /usrx/local/grads/bin/2.0.a3/xgrads -blc "run ${SENSITIVITY}/sen_eof${inum}.gs $center $ndate $inum $valid $lon1 $lon2 $lat1 $lat2 "
  rm cbarn.gs
  rm color.gs
  rm ${SENSITIVITY}/time${inum}norm.ctl
  rm ${SENSITIVITY}/sen_eof${inum}.gs
  rm *.dat
endif

end                                                             ### End loop 4 for three PCs


cd $OUTPUT/
mkdir $ndate
cd $ndate
foreach inum (1 2 3 )                                           ### Begin loop 5 for three z500 PCs
cp ${JOBHOME}cbarn.gs   $OUTPUT/${ndate}/cbarn.gs

if ( $e2 == "21" ) then                                         ###If for colorbar, only 95% significant value will be plot 
cp ${JOBHOME}color.gs  $OUTPUT/${ndate}/color.gs
   else
if ($e2 == "51" ) then
cp ${JOBHOME}color-ec.gs  $OUTPUT/${ndate}/color.gs
endif
if ($e2 == "40") then
cp  ${JOBHOME}color-all.gs  $OUTPUT/${ndate}/color.gs
endif  
endif                                                           ###Endif for colorbar

if ( $e2 == "21" | $e2 == "51" ) then
sed s/__VARI/"${GET_EOF_sed}"time${inum}norm.dat/g $JOBHOME"time_template.ctl" > $SENSITIVITY/time${inum}norm.ctl
echo "XDEF 1 linear 275.0 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
echo "YDEF 1 linear 30 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
cat ${JOBHOME}/member.ctl >> $SENSITIVITY/time${inum}norm.ctl

cp $SENSITIVITY/sen_template_z500.gs  $SENSITIVITY/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/113/${e2}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/999/${ttt}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/timenorm.ctl/"${SENSITIVITY_sed}"time${inum}norm.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/z500_all.ctl/"${GET_BIN_sed}"z500_all_${center}.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${GET_BIN}/z500_all_${center}_samp.ctl | sed -e "s/z500_all_${center}.dat/"${GET_BIN_sed}"z500_all_${center}.dat/g" >file.temp && mv file.temp ${GET_BIN}/z500_all_${center}.ctl
/usrx/local/grads/bin/2.0.a3/xgrads -blc "run ${SENSITIVITY}/sen_eof${inum}_z500.gs $center $ndate $inum $valid $lon1 $lon2 $lat1 $lat2"

rm cbarn.gs
rm color.gs
rm ${SENSITIVITY}/time${inum}norm.ctl
rm ${SENSITIVITY}/sen_eof${inum}_z500.gs
rm *.dat
else
sed s/__VARI/"${GET_EOF_sed}"time${inum}norm.dat/g $JOBHOME"time_template.ctl" > $SENSITIVITY/time${inum}norm.ctl
echo "XDEF 1 linear 275.0 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
echo "YDEF 1 linear 30 ${reso}" >> $SENSITIVITY/time${inum}norm.ctl
cat ${JOBHOME}/member.ctl >> $SENSITIVITY/time${inum}norm.ctl

cp $SENSITIVITY/sen_template_z500_all.gs  $SENSITIVITY/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/113/${e2}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/999/${ttt}/g" > file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/timenorm.ctl/"${SENSITIVITY_sed}"time${inum}norm.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/z500_all.ctl/"${GET_BIN_sed}"z500_all_${center}.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${SENSITIVITY}/sen_eof${inum}_z500.gs | sed -e "s/z500_ano.ctl/"${GET_BIN_sed}"z500_all_${center}_ano.ctl/g" >file.temp && mv file.temp ${SENSITIVITY}/sen_eof${inum}_z500.gs
cat ${GET_BIN}/z500_all_${center}_ano_samp.ctl | sed -e "s/z500_all_${center}_ano.dat/"${GET_BIN_sed}"z500_all_${center}_ano.dat/g" >file.temp && mv file.temp ${GET_BIN}/z500_all_${center}_ano.ctl
/usrx/local/grads/bin/2.0.a3/xgrads -blc "run ${SENSITIVITY}/sen_eof${inum}_z500.gs $center $ndate $inum $valid $lon1 $lon2 $lat1 $lat2"
rm cbarn.gs
rm color.gs
rm ${SENSITIVITY}/time${inum}norm.ctl
rm ${SENSITIVITY}/sen_eof${inum}_z500.gs
rm *.dat
endif
end                                                             ### End loop 5 for three z500 PCs

mv SEN*.gif  ${OUTPUT}/${ndate}/"region"${region}/"day"${day}/
rm color-var.gs
end                                                             ### End loop 3 for valid day
end                                                             ### End loop 2 for region 
end                                                             ### End loop 1 for center
######################################################################
#######If data not available, try 10 times after each 30minutes#######
######################################################################
cd $OUTPUT/
 scp -pr $ndate wd20jd@emcrzdm.ncep.noaa.gov:/home/people/emc/ftp/mmb/sref/CSTAR
    set nn=11
    date
    exit
  else
    @ nn = $nn + 1
    sleep 1800                 ###Yan Comment: this system must use second; 30minutes=30x60seconds.
                   
  endif
echo $nn
end
date
exit

