#!/bin/tcsh
#
# Run quasi-isentropic back-trajectory (QIBT) code
# set echo
# set verbose
set prid = $$
set gribmap = gribmap
set grads = grads 
#
# Set up month and julian day arrays
#
set months = ( jan feb mar apr may jun jul aug sep oct nov dec )
set monls = ( 31 28 31 30 31 30 31 31 30 31 30 31 )
set julians = ( 31 59 90 120 151 181 212 243 273 304 334 365 )
set mons = ( 01 02 03 04 05 06 07 08 09 10 11 12 )
#
# wcr.in contains the following information (as space-delimited ASCII text)
#
# 1 yr      (I)   = Year under consideration (1979 <= yr <= 2003)
# 2 pentst  (I)   = Starting pentad for this job (1 <= pentst <= pentend)
# 3 pentend (I)   = Final pentad for this job (pentst <= pentend <= 73)
# 4 pentad  (I/O) = Current pentad under consideration (pentst <= pentad <= pentend)
# 5 slack   (I)   = Maximum number of days outside (before) pentad to track parcels
# 6 dtmin   (I)   = Trajectory integration time step (minutes)
# 7 juls    (O)   = Julian day at beginning of trajectory window 
# 8 leap    (O)   = Is this pentad after a leap day in a leap year (0/1)?
# 9 leapday (O)   = Does this pentad contain a leap day (0/1)?
#10 pqmin   (I)   = Minimum threshold for parcel water content (%)
#11 partot  (I)   = Number of parcels to launch in pentad
#12 maxcpu  (I)   = Number of parallel processes to run at one time
#13 mstat   (I/O) = Job status flag (set to 0 to start, 1=running or crashed, 2=completed successfully)
#
# (I) must be set before running script; (O) script will set or change this variable
#
set f = `cat wcr.in`
set yr = $f[1]
set pentst = $f[2]
set pentend = $f[3]
set pentad = $f[4]
set slack = $f[5]
set dtmin = $f[6]
set pqmin = $f[10]
set partot = $f[11]
set maxcpu = $f[12]
set mstat = $f[13]
#
# Path and filenames of input regridded precipitation files
#
set stdqcpath = IndianPrecip.regridded
set stdqc = ${stdqcpath}/india_zeroed.${yr}.grd
#
# Output path 
#
set outpath = output.${yr}
#
# Link stdqc (precip) input files
#
rm -f stdqc/stdqc_trp.grd
if ( -e ${stdqc} ) then
  echo stdqc = ${stdqc}
else
  echo ERROR: ${stdqc} not found. Check path and filename
  exit 1
endif
#
cd stdqc
cp ../${stdqc} stdqc_trp.grd
cd ..
#
# Create output path
#
echo output = ${outpath}
mkdir ${outpath}
rm -f output
ln -s ${outpath} output
#
# Avoid running after a job crash
#
 if ( $mstat != 0 ) exit $mstat
#
# Write out wcr.in to show job is in running mode
#
set mstat = 1
cat << -- >! wcr.in
$yr $pentst $pentend $pentad $slack $dtmin $f[7] $f[8] $f[9] $pqmin $partot $maxcpu $mstat
--
#
# Leap year check
#
set leapst = 0
set leapdayst = 0
set leapend = 0
set leapdayend = 0
@ iof = $yr % 4
if ( $yr % 100 == 0 & $yr % 400 != 0 ) set iof = 4
if ( $iof == 0 ) set julians = ( 31 60 91 121 152 182 213 244 274 305 335 366 )
if ( $pentst >= 12 & $iof == 0 ) @ leapst++
if ( $pentst == 12 & $iof == 0 ) @ leapdayst++
if ( $pentend >= 12 & $iof == 0 ) @ leapend++
if ( $pentend == 12 & $iof == 0 ) @ leapdayend++
#
#
echo Proc $prid -- Starting Pentad Loop 
#
##### Pentad Loop ############################################3
#
while ($pentad <= $pentend)
  set leap = 0
  set leapday = 0
  if ( $pentad >= 12 & $iof == 0 ) @ leap++
  if ( $pentad == 12 & $iof == 0 ) @ leapday++
#
  @ yearslack = $yr
  @ yearend = $yr
  @ julend = $pentad * 5 + $leap
  @ julstart = $julend - ( 4 + $leapday )
  set mm = 1
  while ($mm)
    if ($julend <= $julians[$mm]) break
    @ mm++
  end
  @ monend = $mm
  @ julslack = $julstart - $slack
  if ( $julslack < 1 ) then
    @ yearslack = $yearslack - 1
    @ julslack = $julslack + $julians[12]
  endif
  set mm = 1
  set julm1 = 0
  while ($mm)
    if ($julslack <= $julians[$mm]) break
    @ julm1 = $julians[$mm]
    @ mm++
  end
#
# Update wcr.in with leap day information
#
cat << -- >! wcr.in
$yr $pentst $pentend $pentad $slack $dtmin $julslack $leap $leapday $pqmin $partot $maxcpu $mstat
--
#
  @ monslack = $mm
#
# Make appropriate GrADS .ctl files and gribmaps
# to navigate just the data we need
#
  if ($yearslack > 2000) then
    set sgb = 'sig'
  else
    set sgb = 'sgb'
  endif
#
  @ ntsteps = ( 5 + $slack + $leapday ) * 4 + 1
  @ fday = $julslack - $julm1
  set month = $months[$monslack]
  echo Input File Path: flx.${yearslack}
  sed -e "s/%year/${yearslack}/" < flx-ft06.ctl.template >! temp
  sed -e "s/%month/${month}/" < temp >! temp2
  sed -e "s/%ntsteps/${ntsteps}/" < temp2 >! temp
  sed -e "s/%day/${fday}/" < temp >! flx-ft06.ctl
  /bin/rm temp*
#
  if ( -d flx.${yearslack} ) then
    echo Linking flx.${yearslack}
#    ${gribmap} -0 -i flx-ft06.ctl
  else
    echo ERROR: flx.${yearslack} not found.
    exit 1
  endif
#
  echo Input File Path: ${sgb}.${yearslack}
  sed -e "s/%year/${yearslack}/" < sgb-anl.ctl.template >! temp
  sed -e "s/%month/${month}/" < temp >! temp2
  sed -e "s/%ntsteps/${ntsteps}/" < temp2 >! temp
  sed -e "s/%sgb/${sgb}/g" < temp >! temp2
  sed -e "s/%day/${fday}/" < temp2 >! sgb-anl.ctl
  /bin/rm temp*
#
  if ( -d ${sgb}.${yearslack} ) then 
    echo Linking ${sgb}.${yearslack}
#    ${gribmap} -0 -i sgb-anl.ctl
  else
    echo ERROR: ${sgb}.${yearslack} not found.
    exit 1
  endif
#
# Use GrADS to subset GRIB data and write direct-access binaries
# to simplify access to the reanalysis data by the FORTRAN code
#
  ${grads} -blc "run subset_rean.gs"
#
# Run The QIBT code
# Kill the job if the code crashes
#
  echo Running QIBT code for pentad $pentad
  ./qibt2 >! qibt.${pentad}.out || kill -9 $prid
  @ pentad++
#
end
#
##### End Pentad Loop #########################################
#
@ pentad--
set mstat = 2
#
# Successful completion of the job
#
cat << -- >! wcr.in
$yr $pentst $pentend $pentad $slack $dtmin $julslack $leap $leapday $pqmin $partot $maxcpu $mstat
--
#