#!/bin/sh
set -x
#set -nv

cd $DATA

set +x
#########################################################################
echo "----------------------------------------------------------------"
echo "hysplit_conc.sh - Script To Run Hysplit Concentrations and plot."
echo "----------------------------------------------------------------"

# REVISED 19 May 2005 - hourly output cump2 for volcanoes
#         17 Oct 2005 - number particles to MAPTEXT.CFG; MAXPAR
#         19 Jun 2006 - 2nd conc grid for RSMC
#         01 Dec 2006 - zoom and relcen veriables
#         31 Jan 2008 - Google Earth
#         01 Feb 2008 - set end sample date/time to zeros
#                               source attribution (backward)
#         28 Mar 2008 - lower MAXPAR; change filename cdump1 to cdump
#         01 Apr 2008 - change issued to created
#         09 Apr 2008 - force units Bq in MAPTEXT for RSMC 
#         21 Oct 2009 - GDAS-GFS labeling
#         27 Apr 2010 - fix month label in output plot (GDAS-GFS labeling)
#         04 Jun 2010 - multi-processor version
#                       parhplot, effective area source
#         26 May 2011 - error checking
#########################################################################
set -x

##################################
# typenum=1 RSMC
# typenum=2 Volcano
# typenum=3 Homeland Security
# typenum=4 Back-tracking (Source Attribution)
##################################

#############################################################
# Set Output Name for Dispersion/Concentration Simulation
#############################################################
#NOTE FOR HS run:
#       There will be 2 postscript files created (concplot1 and concplot2)
#       concplot1 will contain the fine grid concentration plots, while
#       concplot2 will contain the coarse grid concentration plots

#       If the GIS option (-a) is turned on (1) for concplot then output 
#       GIS text files will also be created.  Each concentration level will 
#       have the following format....
#
#         GIS_level1-level2_(fine or coarse)_(count).txt
#
#       These files can be converted to shapefiles by ascii2shp and input to GIS systems.

#NOTE FOR volcanic ash run:
#       There will be 9 postscript files created 
#         volcplot06.r0 (ash reduction - none)		 6-hour output
#         volcplot06.r1 (ash reduction - small)
#         volcplot06.r2 (ash reduction - medium)
#         volcplot06.r3 (ash reduction - large)
#         volcplot12.r0                       	 	12-hour output for wafs
#         volcplot12.r1 
#         volcplot12.r2
#         volcplot12.r3
#         volcplotae.r0 (alert-ended, 12-hour output for wafs, .rN is required)

#NOTE FOR rsmc run:
#       There will be 1 postscript file created  (concplot.ps)
#############################################################

# run duration above which should run 3-D particle model
# --automatically set when set run duration, and checked at end of setup_hysplit_script
# --longrun and neglongrun must be identical in setup_hysplit_script
# --  and hysplit_conc.sh
# --test is: if [ ${nhrs} -ge ${longrun} -o ${nhrs} -le ${neglongrun} ]
longrun=84 
neglongrun=`echo "${longrun}*-1" | bc`

   ############################################################
   # Set Up Control File for each type of run
   ############################################################

# common grid specs
   ngrd=1
   if [ $dhr2 -ne 0 ]; then ngrd=2; fi

 # grid #1
   gctr[1]="$cnlat $cnlon" 	# center
   gspc[1]="$dlat $dlon"	# spacing
   gspn[1]="$splat $splon"	# span
   gdir[1]="$DATA/"		# directory
 # expand grid for long-range runs
   if [ ${nhrs} -ge ${longrun} -o ${nhrs} -le ${neglongrun} ]; then
      gspn[1]="181. 360."	# span
   fi

 # grid #2
   gctr[2]="$cnlat2 $cnlon2"	# center
   gspc[2]="$dlat2 $dlon2"	# spacing
   gspn[2]="$splat2 $splon2"	# span
   gdir[2]="$DATA/"		# directory
   gfil[2]="cdump2"		# filename

case $typenum in 
1) # RSMC - particle with wet and dry deposition and half-life
   #        1 concentration grid
   npol=1
   pol[1]=$ident
   qrate[1]=$qrat

	# grid specs
   gfil[1]="cdump"		# filename
   gnlv[1]=2	 		# number of levels
   glvl[1]="0 $height"		# levels

   year=`echo $csyr | cut -c3-`
 # RSMC protocol - samples start at 00z or 12z
   if [ $dshr -ge 12 ]
   then
      gsst[1]="$year $csmo $csda 12 00"	
   else
      gsst[1]="$year $csmo $csda 00 00"
   fi
   gsnd[1]="00 00 00 00 00"	# sample stop
   gint[1]="00 $dhr 00"		# sampling interval 

   if [ $ngrd -eq 2 ]; then

 # number of hours of output  (from present time)
   dfcst=18

 # round rlse down to previous multiple of $dhr2 as start time (3 for RSMC)
     # release start time
       YR=$dsyr 
       MO=$dsmo
       DA=$dsda
       HR=$dshr
     # sample start time = release start time rounded down
       if [ ! `expr $HR % $dhr2` -eq 0 ]; then
          HR=`expr $HR - $HR % $dhr2`
          if [ $HR -lt 10 ]; then HR="0"$HR; fi
       fi 

     # present time
       PYR=`date -u +%Y`
       PMO=`date -u +%m`
       PDA=`date -u +%d`
       PHR=`date -u +%H`
     # sample end time = (present + 18 hours) rounded up
       if [ `expr $PHR % $dhr2` -eq 0 ]; then
          inc=$dfcst
          NS=`expr $dfcst / $dhr2` 		# number of samples
       else
          inc=`expr $PHR % $dhr2`
          inc=`expr $dfcst + $dhr2 - $inc`
          NS=`expr $dfcst / $dhr2 + 1`
       fi 
       END=`${utilexec}/ndate $inc ${PYR}${PMO}${PDA}${PHR}`
       EYR=`echo $END | cut -c1-4`
       EMO=`echo $END | cut -c5-6`
       EDA=`echo $END | cut -c7-8`
       EHR=`echo $END | cut -c9-10`

     # time periods (samples) N1 to N2 will be plotted
       NHR=`${utilexec}/nhour ${EYR}${EMO}${EDA}${EHR} ${YR}${MO}${DA}${HR}`
       N2=`expr $NHR / $dhr2`
       N1=`expr $N2 - $NS + 1`

########################################################################
   gnlv[2]=2			# number of levels
   glvl[2]="0 $height2"		# levels
   gsst[2]="$YR $MO $DA $HR 00"		# sample start 
   gsnd[2]="$EYR $EMO $EDA $EHR 00"	# sample stop 
   gint[2]="00 $dhr2 00"	# sampling interval 	

   fi

	# depositing specs
   pspc[1]="1.0 1.0 1.0"		# particle specs
   vdep[1]="$dryvl 0.0 0.0 0.0 0.0"	# deposition velocity etc.
   wet[1]="0.0 $wetin $wetlo"		# wet removal
   half[1]=$rhalf 			# half-life
   res[1]="0.0"				# resuspension
   ;;

2) # Volcano - 4 particle sizes, with gravitational settling
   #           2 concentration grids; 1 6-hourly, 2nd 1-hourly
   #           sample start time and run time depend on eruption time
   #             with respect to meteo initialization time
   npol=4
   if [ ${qrat} -eq 0 ]
   then
      pol[1]=p006; qrate[1]=0.0   		#  0.3 -  1 microns
      pol[2]=p020; qrate[2]=0.0   		#  1   -  3 microns
      pol[3]=p060; qrate[3]=0.0   		#  3   - 10 microns
      pol[4]=p200; qrate[4]=0.0   		# 10   - 30 microns
   else
      pol[1]=p006; qrate[1]=0.008 		#  0.3 -  1 microns
      pol[2]=p020; qrate[2]=0.068 		#  1   -  3 microns
      pol[3]=p060; qrate[3]=0.25  		#  3   - 10 microns
      pol[4]=p200; qrate[4]=0.67  		# 10   - 30 microns
   fi

 # total emission is 1.0
   if [ ${qrat} -ne 0 ]
   then
      for i in 1 2 3 4
      do
        qrate[${i}]=`echo "scale=5\n ${qrate[${i}]}/${qhrs}" | bc`
      done
   fi

	# grid specs for conc grid #1
   gfil[1]="cdump"		# filename
   gnlv[1]=3			# number of levels
   glvl[1]="6096 10668 16764"	# levels (needed for volcplot)
   gsst[1]="00 00 00 00 00"	# sample start specified below
   gsnd[1]="00 00 00 00 00"	# sample stop	
   gint[1]="-1 6 00"	# sampling interval  - hourly concentration every 6-h


 # start time is met fcst initialization time
  #YR=`echo ${PDY} | cut -c1-4`
  #MO=`echo ${PDY} | cut -c5-6`
  #DA=`echo ${PDY} | cut -c7-8`
  #HR=${cyc}

# grid specs for conc grid #2 (for SAB NAWIPS)

 # start time is the time of the eruption, if up to 6 hours ago,
 #        or, is 6 hours ago, if the eruption was more than 6 hours ago
  
 # present time
   YR=`date -u +%Y`
   MO=`date -u +%m`
   DA=`date -u +%d`
   HR=`date -u +%H`
   
 # 6 hours before present
   back=6
   inc=`expr -1 \* $back` 
   SIXBP=`${utilexec}/ndate $inc ${YR}${MO}${DA}${HR}`
   syr=`echo $SIXBP | cut -c1-4`
   smo=`echo $SIXBP | cut -c5-6`
   sda=`echo $SIXBP | cut -c7-8`
   shr=`echo $SIXBP | cut -c9-10`

 # eruption start time
    # $dsyr 
    # $dsmo
    # $dsda
    # $dshr

 # desired number of forecast hours for hourly output
   dfcst=18

 # extra hours to have if needed 
   extra=3

 # dt = difference in hours between present and eruption start
 # dt>0 ==> eruption before present
     dt=`${utilexec}/nhour ${YR}${MO}${DA}${HR} $dsyr$dsmo$dsda$dshr`
     let dt=${dt}		# to get rid of leading "0"
   
 # dhour = duration of cdump2 sampling 
   if [ ${dt} -ge ${back} ]; then
    # eruption starts at least X hours before present(BP), start hourlies 6h BP
    # grid #2 sample start time
      YR=$syr
      MO=$smo
      DA=$sda
      HR=$shr
    # grid #2 sample end time
      let dhour=${dfcst}+${extra}+${back}
   elif [ ${dt} -ge 0 ]; then
    # eruption started at the present or up to 6h BP, start hourlies at eruption
    # grid #2 sample start time
      YR=$dsyr
      MO=$dsmo
      DA=$dsda
      HR=$dshr
    # grid #2 sample end time
      let dhour=${dfcst}+${extra}+${dt}
   else     
    # don't do hourly plots for eruption starting after present (not real erupt)
  #   ngrd=1
    # OR set start time to eruption time and output hourlies for 18 hours only
    # grid #2 sample start time
      YR=$dsyr
      MO=$dsmo
      DA=$dsda
      HR=$dshr
    # grid #2 sample end time
      dhour=${dfcst}
   fi

 # grid #2 sample end time
   echo $dhour >> DHOUR
   SMPET=`${utilexec}/ndate $dhour ${YR}${MO}${DA}${HR}`
   EYR=`echo $SMPET | cut -c1-4`
   EMO=`echo $SMPET | cut -c5-6`
   EDA=`echo $SMPET | cut -c7-8`
   EHR=`echo $SMPET | cut -c9-10`

   gnlv[2]=3			# number of levels
   glvl[2]="6096 10668 16764"	# levels (want same levels in NAWIPS)
   gsst[2]="$YR $MO $DA $HR 00"		# sample start 
   gsnd[2]="$EYR $EMO $EDA $EHR 00"	# sample stop 
   gint[2]="-1 $dhr2 00"	# sampling interval 	

     ###########################################
     # sample start time and run time - volcano
     ###########################################

 # difference between eruption time (Te) and met fcst init time (Tm) in hours
     # DT was calculated in hysplit_trajconc.sh
     # positive for eruption after Tm
     # negative              before

                       ## don't want output plots during analysis period
   # if Te is at or after 3 hours before Tm
   if test $DT -ge -3 
   then
    # sample start time = eruption time
      syr=`echo $dsyr | cut -c3-4`
      smo=$dsmo 
      sda=$dsda
      shr=$dshr

   else

    # if Te is at or after 18-h before Tm and before 3-h before Tm
      if test $DT -ge -18 -a $DT -lt -3
      then
       # sample start time = Te + 12 
         inc=12
         SMPST=`${utilexec}/ndate $inc $dsyr$dsmo$dsda$dshr`
         syr=`echo $SMPST | cut -c3-4`
         smo=`echo $SMPST | cut -c5-6`
         sda=`echo $SMPST | cut -c7-8`
         shr=`echo $SMPST | cut -c9-10`

      elif [ $DT -lt -174 ]; then
      #  could handle earlier eruptions by increasing values for N 
      #   as long as have meteorology to support it
         msg='ERROR This run is starting too far back in time.'
         postmsg "$jlogfile" "$msg"
         export err=1; err_chk 

      else
        for N in 3 4 5 6 7 8 9 10 11 12 13 14 15
        do
            # for next upgrade - consider simplifying this by just
            # setting syr/smo/sda/shr to the current time

                                                # N=3: -30 --> -18 hours
                                                # N=4: -42 --> -30 hours
                                                # N=5: -54 --> -42 hours
                                                # N=6: -66 --> -54 hours
                                                # N=7: -78 --> -66 hours
                                                # N=8: -90 --> -78 hours
                                                # N=9: -102 --> -90 hours
                                                # N=10:-114 -->-102 hours
                                                # N=11:-126 -->-114 hours
                                                # N=12:-138 -->-126 hours
                                                # N=13:-150 -->-138 hours
                                                # N=14:-162 -->-150 hours
                                                # N=15:-174 -->-162 hours

           lo=`echo "-12*$N+6" | bc`
           hi=`echo "-12*($N-1)+6" | bc`
  
           if test $DT -ge $lo -a $DT -lt $hi
           then
             inc=`echo "12*($N-1)"	| bc`
             SMPST=`${utilexec}/ndate $inc $dsyr$dsmo$dsda$dshr`
             syr=`echo $SMPST | cut -c3-4`
             smo=`echo $SMPST | cut -c5-6`
             sda=`echo $SMPST | cut -c7-8`
             shr=`echo $SMPST | cut -c9-10`
           fi
        done

      fi
   fi

   gsst[1]="$syr $smo $sda $shr 00"	# sample start

	# depositing specs
   pspc[1]="0.6 2.5 1.0"	# particle specs (size in microns, dens, shape)
   vdep[1]="0.0 0.0 0.0 0.0 0.0"	# deposition velocity etc.
   wet[1]="0.0 ${wetin} ${wetlo}"	# wet removal
   half[1]="0.0"		# half-life
   res[1]="0.0"			# resuspension

   pspc[2]="2.0 2.5 1.0"	
   vdep[2]=${vdep[1]}     
   wet[2]=${wet[1]}	
   half[2]=${half[1]}
   res[2]=${res[1]} 

   pspc[3]="6.0 2.5 1.0"	
   vdep[3]=${vdep[1]}     
   wet[3]=${wet[1]}	
   half[3]=${half[1]}
   res[3]=${res[1]} 

   pspc[4]="20.0 2.5 1.0"	
   vdep[4]=${vdep[1]}     
   wet[4]=${wet[1]}	
   half[4]=${half[1]}
   res[4]=${res[1]} 
   ;;

3) # Homeland Security - gaseous pollutant, non-depositing
   #                     2 concentration grids: 1st 12-hours, after 12-h
   npol=1
   pol[1]=$ident
   qrate[1]=$qrat
                           # ------------------------------------
                           # fine grid output 1st $delhr hours
                           # ------------------------------------
	# grid specs
   gfil[1]="cdump"		# filename
   gnlv[1]=2			# number of levels
   glvl[1]="10 $height"		# levels
   year=`echo $csyr | cut -c3-`
   gsst[1]="$year $csmo $csda $cshr 00"	# sample start
   gsnd[1]="00 00 00 $delhr 00"	# sample end
   gint[1]="00 $dhr 00"		# sampling interval	     
                           # -----------------------------------------
                           # coarse grid output beyond $delhr hours
                           # -----------------------------------------
   if [ $ngrd -eq 2 ]; then
   gnlv[2]=2			# number of levels
   glvl[2]="10 $height2"	# levels
   gsst[2]="00 00 00 $delhr 00"	# sample start
   gsnd[2]="00 00 00 00 00"	# sample end
   gint[2]="00 $dhr2 00"	# sampling interval 	
   fi

	# depositing specs
   if [ ${wetin} -eq 0 -a ${wetlo} -eq 0 ]
   then
      pspc[1]="0.0 0.0 0.0"		# particle specs  0's --> gas
   else
      pspc[1]="1.0 1.0 1.0"		# particle specs  1's --> particle
   fi
   vdep[1]="$dryvl 0.0 0.0 0.0 0.0"	# deposition velocity etc.
   wet[1]="0.0 $wetin $wetlo"		# wet removal
   half[1]="$rhalf"			# half-life
   res[1]="0.0"				# resuspension
   ;;
4) # Back-tracking (Source Attribution)
   #        no deposition, 1 concentration grid, gas (not particle)
   npol=1
   pol[1]=$ident
   qrate[1]=$qrat

	# grid specs
   gfil[1]="cdump"		# filename
   gnlv[1]=1	 		# number of levels
   glvl[1]="$height"		# levels

   year=`echo $csyr | cut -c3-`
   gsst[1]="$year $csmo $csda $cshr 00"
   gsnd[1]="00 00 00 00 00"	# sample stop
   gint[1]="00 $dhr 00"		# sampling interval 

 # number of hours of output  (from present time)
   dfcst=18

 # round rlse down to previous multiple of $dhr2 as start time (3 for RSMC)
     # release start time
       YR=$dsyr 
       MO=$dsmo
       DA=$dsda
       HR=$dshr
     # sample start time = release start time rounded down
      #if [ ! `expr $HR % $dhr2` -eq 0 ]; then
      #   HR=`expr $HR - $HR % $dhr2`
      #   if [ $HR -lt 10 ]; then HR="0"$HR; fi
      #fi 

     # present time
       PYR=`date -u +%Y`
       PMO=`date -u +%m`
       PDA=`date -u +%d`
       PHR=`date -u +%H`
     # sample end time = (present + 18 hours) rounded up
      #if [ `expr $PHR % $dhr2` -eq 0 ]; then
      #   inc=$dfcst
      #   NS=`expr $dfcst / $dhr2` 		# number of samples
      #else
      #   inc=`expr $PHR % $dhr2`
      #   inc=`expr $dfcst + $dhr2 - $inc`
      #   NS=`expr $dfcst / $dhr2 + 1`
      #fi 
     # END=`${utilexec}/ndate $inc ${PYR}${PMO}${PDA}${PHR}`
     # EYR=`echo $END | cut -c1-4`
     # EMO=`echo $END | cut -c5-6`
     # EDA=`echo $END | cut -c7-8`
     # EHR=`echo $END | cut -c9-10`
     # force to zeros for now - need to calc based on request/avail met
       EYR=00
       EMO=00
       EDA=00
       EHR=00

     # time periods (samples) N1 to N2 will be plotted
      #NHR=`${utilexec}/nhour ${EYR}${EMO}${EDA}${EHR} ${YR}${MO}${DA}${HR}`
      #N2=`expr $NHR / $dhr2`
      #N1=`expr $N2 - $NS + 1`

########################################################################
   gnlv[2]=2			# number of levels
   glvl[2]="0 $height2"		# levels
   gsst[2]="$YR $MO $DA $HR 00"		# sample start 
   gsnd[2]="$EYR $EMO $EDA $EHR 00"	# sample stop 
   gint[2]="00 $dhr2 00"	# sampling interval 	

	# depositing specs
   pspc[1]="0.0 0.0 0.0"		# particle specs
   vdep[1]="$dryvl 0.0 0.0 0.0 0.0"	# deposition velocity etc.
   wet[1]="0.0 $wetin $wetlo"		# wet removal
   half[1]=$rhalf 			# half-life
   res[1]="0.0"				# resuspension

   ;;

esac

  ########################
  # write CONTROL file (this is what hysplit reads)
  ########################

  echo "$ibyr $ibmo $ibda $ibhr" >CONTROL.conc	# calculation start
  if [ ${olat2} -eq -99 -o ${olon2} -eq -99 -o ${olat3} -eq -99 -o ${olon3} -eq -99 ]
  then
     # this test is also in hysplit_post.sh make open symbol on output plot
  echo "2                      ">>CONTROL.conc	# number of starting locations
  echo "$olat $olon ${olvl1%.*}">>CONTROL.conc	# loc #1 lat lon height
  echo "$olat $olon ${olvl2%.*}">>CONTROL.conc	# loc #2 lat lon height
  else
  echo "6                      ">>CONTROL.conc	# number of starting locations
  echo "$olat $olon ${olvl1%.*}">>CONTROL.conc	# loc #1 lat lon height
  echo "$olat $olon ${olvl2%.*}">>CONTROL.conc	# loc #2 lat lon height
  echo "$olat2 $olon2 ${olvl1%.*}">>CONTROL.conc	# loc #1 lat lon height
  echo "$olat2 $olon2 ${olvl2%.*}">>CONTROL.conc	# loc #2 lat lon height
  echo "$olat3 $olon3 ${olvl1%.*}">>CONTROL.conc	# loc #1 lat lon height
  echo "$olat3 $olon3 ${olvl2%.*}">>CONTROL.conc	# loc #2 lat lon height
  fi
  echo "$nhrs                  ">>CONTROL.conc	# total run time
  echo "0                      ">>CONTROL.conc	# vertical motion 
  echo "$zdata	               ">>CONTROL.conc	# top of model domain
  echo "$nfile"                 >>CONTROL.conc	# number of input meteo files
                             # meteo files may appear in any order
                             # hysplit will choose appropriate files in order
                               # see MESSAGE file
                             # meteo files must exist
  ifile=0
  if [ ! "$met1" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/                 ">>CONTROL.conc	# file #1 directory
     echo "${met1}              ">>CONTROL.conc	# file #1 filename
  fi
  if [ ! "$met2" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #2 directory
     echo "${met2}           ">>CONTROL.conc	# file #2 filename
  fi
  if [ ! "$met3" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #3 directory
     echo "${met3}           ">>CONTROL.conc	# file #3 filename
  fi
  if [ ! "$met4" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #4 directory
     echo "${met4}           ">>CONTROL.conc	# file #4 filename
  fi
  if [ ! "$met5" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #5 directory
     echo "${met5}           ">>CONTROL.conc	# file #5 filename
  fi
  if [ ! "$met6" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #5 directory
     echo "${met6}           ">>CONTROL.conc	# file #5 filename
  fi
  if [ ! "$met7" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #5 directory
     echo "${met7}           ">>CONTROL.conc	# file #5 filename
  fi
  if [ ! "$met8" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #5 directory
     echo "${met8}           ">>CONTROL.conc	# file #5 filename
  fi
  if [ ! "$met9" = "undefined" ]
  then
     let ifile=$ifile+1
     echo "$DATA/              ">>CONTROL.conc	# file #5 directory
     echo "${met9}           ">>CONTROL.conc	# file #5 filename
  fi

  if [ $ifile -ne $nfile ]
  then
     msg='ERROR - Unexpected number of meteorology files written to CONTROL.conc - hymodelm will fail'
     postmsg "$jlogfile" "$msg"
     export err=2; err_exit 
  fi
  echo "$npol                  ">>CONTROL.conc	# number of pollutants
  i=1
  while [ $i -le $npol ]; do
    echo "${pol[$i]}                ">>CONTROL.conc   # pollutant ID 
    echo "${qrate[$i]}              ">>CONTROL.conc   # emission rate (per hour)
    echo "$qhrs                     ">>CONTROL.conc   # hours of emission
    year=`echo $dsyr | cut -c3-`
    echo "$year $dsmo $dsda $dshr $dsmn">>CONTROL.conc   # release start 
    let i=i+1
  done
  echo "$ngrd           ">>CONTROL.conc	# number of simultaneous conc grids
  i=1
  while [ $i -le $ngrd ]; do
    echo "${gctr[$i]}          ">>CONTROL.conc	# grid i center lat lon
    echo "${gspc[$i]}          ">>CONTROL.conc	# grid i grid spacing lat lon 
    echo "${gspn[$i]}          ">>CONTROL.conc	# grid i grid span lat lon
    echo "${gdir[$i]}          ">>CONTROL.conc	# grid i directory
    echo "${gfil[$i]}          ">>CONTROL.conc	# grid i filename
    echo "${gnlv[$i]}          ">>CONTROL.conc	# grid i number of levels
    echo "${glvl[$i]}          ">>CONTROL.conc	# grid i levels
    echo "${gsst[$i]}          ">>CONTROL.conc	# grid i sampling start
    echo "${gsnd[$i]}          ">>CONTROL.conc	# grid i sampling end 
    echo "${gint[$i]}          ">>CONTROL.conc	# grid i sampling interval
    let i=i+1
  done
  echo "$npol                  ">>CONTROL.conc	# number of pollutants depositing

  i=1
  while [ $i -le $npol ]; do
    echo "${pspc[$i]}          ">>CONTROL.conc	# particle specs 
    echo "${vdep[$i]}          ">>CONTROL.conc	# deposition velocity etc
    echo "${wet[$i]}           ">>CONTROL.conc	# wet removal 
    echo "${half[$i]}          ">>CONTROL.conc	# half-life
    echo "${res[$i]}           ">>CONTROL.conc	# resuspension
    let i=i+1
  done

#############################
# now that control file is written, set gfile[2] to cdump for volcplot
#############################
   case $typenum in
     2) gfil[2]=${gfil[1]};;
   esac

###############################################################################
# write particle size specs for volcanic ash
#  for use in hysplit_parshift_setup
###############################################################################
  case $typenum in
     2) echo "# ID   diameter(micrometers, um) " >PARTICLE.SIZES
        i=1
        while [ $i -le $npol ]
        do
        # extract diameter only from $pspc
          echo ${pspc[$i]} | read dia[$i] rest
          echo "${i} ${pol[$i]} ${dia[$i]}">>PARTICLE.SIZES
          let i=${i}+1
        done
        ;;
  esac

###############################################################################
# Set up namelist parameters for each type of Dispersion/Conc Simulation
#   only those different from defaults must be given
###############################################################################

POUTF=PARDUMP

rm -f SETUP.conc
echo "&SETUP"             > SETUP.conc
echo "INITD = $initd,"   >> SETUP.conc	# initial particle distribution
echo "KHMAX = $khmax,"   >> SETUP.conc	# max particle age
echo "NUMPAR = $numpar," >> SETUP.conc	# number of particles
echo "MAXPAR = 50000,"   >> SETUP.conc  # max number of particles
echo "ISOT = $isot,"     >> SETUP.conc	# isotropic turbulence
echo "POUTF = '${POUTF}',">> SETUP.conc	# particle output file
echo "NDUMP = $ndump,"   >> SETUP.conc	# write PARDUMP beginning at this simulation hour
echo "NCYCL = $ncycl,"   >> SETUP.conc	# write PARDUMP at this time interval (hours)

   ############################################################
   # get parinit file from /com if using it (particle position initialization file)
   #   and write filename to SETUP file
   ############################################################
   if [ "${pinpf}" != "" ]
   then
      if [ -s ${pinpf} ]				# either PARDUMP or PARSHIFT
      then
         cp ${pinpf} PARINIT
         echo "PINPF = 'PARINIT',">> SETUP.conc 	# particle input file
      elif [ -s ${pinpf}.001 ]
      then
         cp ${pinpf}.001 PARINIT.001
         cp ${pinpf}.002 PARINIT.002
         cp ${pinpf}.003 PARINIT.003
         cp ${pinpf}.004 PARINIT.004
         cp ${pinpf}.005 PARINIT.005
         cp ${pinpf}.006 PARINIT.006
         cp ${pinpf}.007 PARINIT.007
         cp ${pinpf}.008 PARINIT.008
         echo "PINPF = 'PARINIT',">> SETUP.conc 	# particle input file
      else
         echo "ERROR: Initialization file(s) ${pinpf} not found"
         exit 301
      fi
   fi

echo "KMSL = $kmsl"      >> SETUP.conc  # agl/msl
if [ ${nhrs} -ge ${longrun} -o ${nhrs} -le ${neglongrun} ]; then
   echo "DELT = 20"      >> SETUP.conc  # time step (minutes)
fi 
echo "/"                 >> SETUP.conc

# exit 		# if just want to see CONTROL.conc and SETUP.conc
                # but hysplit_post.sh still runs

############################################################
# list CONTROL and SETUP for debugging
############################################################
set +x
echo "-------------------------"
echo "CONTROL.conc---------------"
set -x
cat CONTROL.conc
set +x
echo "SETUP.conc---------------"
set -x
cat SETUP.conc
set +x
echo "-------------------------"
set -x

exit