#!/bin/bash
#=====================================================================
#
# NAME 
#   xwambdu
#
# PURPOSE
#   Master script to run *ENLIL* with *WSA* daily-updated data.
#
# BACKGROUND
#   ENLIL heliospheric computations are driven by boundary values
#   calculated from the WSA coronal maps as summarized below.
#
#   The coronal maps are placed in the <cordir>/<cordata>/<corobs>
#   directory.  Alternatively, a single coronal map <cordir>/<corfile>
#   can be used. The sub-directory <cordata> = wsa<vers> indicates
#   data from different WSA versions and <corobs> = gongb, gongz.
# 
#   Pre-processing, heliospheric, and post-processing codes are
#   compiled and placed in:
#   - <projdir>/<project>/code.... directory;
#   grid and boundary values (produced by pre-processing codes) are in:
#   - <projdir>/<project>/case.... directory;
#   heliospheric results (produced by enlil and post-processing codes)
#   are in:
#   - <projdir>/<project>/run.... directory.
#
#   The pre-processing "reg2grd" code uses the input files:
#   - reg2grd.in (prepared by this script, with files locations and
#     model-free parameters);
#   and produces the output files:
#   - reg2grd.out (log file);
#   - grd.nc (numerical grid).
#
#   The pre-processing "wsa2bc" code uses the input files:
#   - wsa2bc.in (prepared by this script, with files locations and
#     model-free parameters);
#   - WSA coronal maps (fits files); 
#   and produces the output files:
#   - wsa2bc.out (log file);
#   - bnd.nc (boundary conditions).
#   
#   The heliospheric "enlil" code uses the input files:
#   - enlil.in (prepared by this script, with files locations and
#     model-free parameters);
#   - grd.nc (numerical grid).
#   - bnd.nc (boundary conditions).
#   and produces the output files:
#   - enlil.out (log file);
#   - evo.<namobs>.nc (evolution at planets and spacecraft);
#   - fld.<namobs>.nc (values along IMF lines passing through objects);
#   - tim.<rec>.nc (values in the heliospheric computational domain).
#
#   There are various post-processing codes and IDL visualization procedures
#   as used in this script.
#
# USAGE
#   xwambdu [-<parameter> <value>]
#   xwambdu make [-<parameter> <value>]
#   xwambdu case [-<parameter> <value>]
#   xwambdu plot_case [-<parameter> <value>]
#   xwambdu run [-<parameter> <value>]
#   xwambdu data_run [-<parameter> <value>]
#   xwambdu plot_run [-<parameter> <value>]
#   xwambdu clean [-<parameter> <value>]
#
# DESCRIPTION
#   make      - Set array dimensions and compile the codes
#   case      - Set input parameters and execute the case code
#   plot_case - Plot results of the case computation
#   run       - Set input parameters and execute the run code
#   data_run  - Data results of the run computation
#   plot_run  - Plot results of the run computation
#   clean     - Clean the project directory
#
# PARAMETERS - PROJECT
#   -project <char>   - Project name:
#                       wamb/wambdu
#
# PARAMETERS - NUMERICAL MODEL
#   -numo <char>      - Numerical model setting:
#                       mp1um1d | mp1va2d
#
# PARAMETERS - COMPUTATIONAL REGION AND GRID
#   -nreg <int>       - Computational grid mode:
#                       1 | 2 | 3
#   -reg <int>        - Computational region (AU):
#                       1 | 2 | 3 | 6 | 12
#   -res <char>       - Numerical grid resolution:
#                       low | med | high | fine | super | ultra
#
# PARAMETERS - COMPUTATIONAL PERIOD AND OUTPUT
#   -obsdate <char>   - Use observations at or up to this date:
#                       null | <yyyy-mm-ddThh>
#   -rundate <char>   - Start computation at this date:
#                       2015-08-01T00 | <yyyy-mm-ddThh>
#   -tstop <float>    - Stop computation at this time (h):
#                       168.
#   -tstep <float>    - Output with this step in time (h):
#                       24.
#
# PARAMETERS - BACKGROUND SOLAR WIND
#   -cordata <char>   - Coronal data sub-directory:
#                       wsa22 | wsa45 | wsa52 | wsada | wsafr
#   -corfile <char>   - Coronal data file (null if cordata/corymdh used):
#                       null | <name>
#   -cormode <char>   - Coronal data mode:
#                       single | multi
#   -corobs <char>    - Solar photospheric field observatory:
#                       gongb | gongz | mdi | mwo | nso
#   -corymdh <char>   - Coronal data sub-sub-directory structure by date:
#                       yyyy/mm | yyyy/mm/dd | null
#   -amb <char>       - Ambient wind conditions setting:
#                       a8b1 | a3b2 | a6b1
#   -bndrot <char>    - Rotation of the inner boundary:
#                       synodic | sidereal | null
#   -gamma <float>    - Ratio of specific heats:
#                       1.6666667
#   -xalpha <float>   - Fraction of alpha particles (rel. to protons)
#                       0.
#
# PARAMETERS - RUN
#   -qheat <float>    - Volumetric heating factor
#                       0.3
#
# PARAMETERS - PLOT PROCEDURES
#   -pbnd <int>       - Plot_bnd if non-zero:
#                       1
#   -pbnd2 <int>      - Plot_bnd2 if non-zero:
#                       0
#   -pbnd6 <int>      - Plot_bnd6 if non-zero:
#                       0
#   -pbndtim1 <int>   - Plot_bndtim1 if non-zero:
#                       0
#   -pevo1o5 <int>    - Plot_evo1o5 if non-zero:
#                       0
#   -pevo5o1 <int>    - Plot_evo5o1 if non-zero:
#                       1
#   -ptim1 <int>      - Plot_tim1 with this step in output files:
#                       0
#   -ptim1o5 <int>    - Plot_tim1o5 with this step in output files:
#                       0
#   -ptim2e5 <int>    - Plot_tim2e5 with this step in output files:
#                       0
#   -ptim3 <int>      - Plot_tim3 with this step in output files:
#                       0
#   -ptim3e1 <int>    - Plot_tim3e1 with this step in output files:
#                       1
#
# PARAMETERS - PLOT OPTIONS
#   -fontsize <float> - Font size (1=small, 2=medium, 3=large)
#                       1.0
#   -info <int>       - Informative footnote (1=yes, 0=no, -1=ghost):
#                       1
#   -labitem <char>   - Label plot by item (optional, at top-left):
#                       DEMO
#   -rmax <float>     - Plotting region (AU):
#                       1.7
#   -size <int>       - Size (0=720, 1=960, 2=1440, 3=1920, 4=3840):
#                       1
#   -wide <int>       - Width aspect ratio (0=1:1, 1=1.5:1, 2=2:1 ,3=3:1):
#                       1
#   -zoom <int>       - Zoom (0=none | 1=66% | 2=50% (or 33%):
#                       2
#
#   Note: The first value in the list is the default value.
#
# INSTALLATION
#   You have to edit in the ENLIL distribution:
#   - environment variables in env.inc script.
#   You have to edit in this script:
#   - project name and directory (where results will go);
#   - data directory (repository of coronal maps).
#
# NOTE
#   Default parameters correspond to a low-resolution computation
#   using the reference numerical model and ambient wind settings
#   for the 2015 August period without CMEs in a heliospheric region
#   between 0.1 and 1.1 AU for -168+168 h (-7+7 days).
#
# NOTE
#   It is assumed you are in the <mydir> directory path:
#   .../<HELIO>/src/<ENLIL>/scripts/<mydir>
#   <HELIO_DIR> contains: data, env, fcst, opt, proj, src
#   <ENLIL_DIR> contains: codes, data, scripts, vis
#
#---------------------------------------------------------------------
# SYSTEM DEPENDENT PARAMETERS
#---------------------------------------------------------------------
# Base directories

sname=$(realpath "$0")       # Script full-path name
sdir=$(dirname $sname)       # Directory where the script is located

ENLIL_DIR=${sdir%/*/*/*}     # ENLIL source directory
HELIO_DIR=${sdir%/*/*/*/*/*} # HELIO-WEATHER system directory

#---------------------------------------------------------------------
# Project and data directories

echo source $HELIO_DIR/env/dirs.inc
     source $HELIO_DIR/env/dirs.inc

#---------------------------------------------------------------------
# Operating system

system=linux        # Operating system (linux | aix)
batch=no            # Submit run to batch queue (no | yes)
binsrc=src          # Numerical code (bin | src)
vis=idl             # Visualization software (idl | gdl)

#---------------------------------------------------------------------
# PROJECT DEFAULT PARAMETERS
#---------------------------------------------------------------------
# Default project parameters

project=wamb/wambdu # Project name

#---------------------------------------------------------------------
# Data sub-directory

heldata=helioweb    # Heliospheric data sub-directory

#---------------------------------------------------------------------
# Default code parameters

procs=16            # No. of processors (1 | 2 | 4 | 8 | 16 | 32 | 64)
check=nocheck       # Compile with diagnostic checks (nocheck | check)
outlog=outlog       # Output log file (outlog | noutlog)

#---------------------------------------------------------------------
# Default numerical model parameters

numo=mp1va2d        # Numerical model setting (mp1um1d | mp1va2d)

#---------------------------------------------------------------------
# Default computational region and grid parameters

nreg=1              # Computational grid mode (1 | 2 | 3)
reg=1               # Computational region (AU)
res=low             # Numerical grid resolution
res23=1             # Numerical grid resolution - finer in x2+x3
npos=5              # Number of positions for evolutions (1 | 3 | 5)

#---------------------------------------------------------------------
# Default computational period and output parameters (relative to "rundate")

obsdate=null        # Use observations up to this date (null | yyyy-mm-ddThh)
rundate=2011-11-27  # Start computation at this date (yyyy-mm-ddThh)
tstop=720.          # Stop computation at this time (h)
tstep=720.          # Output with this step in time (h)

#---------------------------------------------------------------------
# Default ambient solar wind parameters

# Input data
cordata=wsa22       # Coronal data sub-directory (wsa22 | wsa45 | wsa52 | wsada | wsafr)
corfile=null        # Coronal data file (null if cordata/corymdh used)
cormode=single      # Coronal data mode (single | multi)
corobs=gongz        # Observatory name (gongb | gongz | mdi | mwo | nso)
corymdh=yyyy/mm     # Coronal data sub-sub-directory structure by date (yyyy/mm | yyyy/mm/dd | null)

# Model parameters
amb=a6b1            # Ambient solar wind conditions setting (a8b1 | a3b2 | a6b1 | ace3a)
bndrot=synodic      # Rotation of the inner boundary (null | sidereal | synodic)
gamma=1.6666667     # Ratio of specific heats
xalpha=0.           # Fraction of alpha particles (rel. to protons)

#---------------------------------------------------------------------
# Default run parameters

cdifb=0.2           # Div(B) diffusion coefficient
qheat=0.3           # Volumetric heating factor

#---------------------------------------------------------------------
# Default plot procedures

pbnd=1              # Plot_bnd if non-zero
pbnd2=0             # Plot_bnd2 if non-zero
pbnd6=0             # Plot_bnd6 if non-zero
pbndtim1=0          # Plot_bndtim1 if non-zero
pevo1o5=0           # Plot_evo1o5 if non-zero
pevo5o1=1           # Plot_evo5o1 if non-zero
ptim1=0             # Plot_tim1 with this step in output files
ptim1o5=0           # Plot_tim1o5 with this step in output files
ptim2e5=0           # Plot_tim2e5 with this step in output files
ptim3=0             # Plot_tim3 with this step in output files
ptim3e1=1           # Plot_tim3e1 with this step in output files

#---------------------------------------------------------------------
# Default plot options

fontsize=1.0        # Font size (1=small, 2=medium, 3=large)
info=1              # Informative footnote (1=yes, 0=no, -1=ghost)
labitem=DEMO        # Label plot by item (optional, at top left)
rmax=1.7            # Plotting region (AU)
size=1              # Size (0=720, 1=960, 2=1440, 3=1920, 4=3840)
wide=1              # Width aspect ratio (0=1:1, 1=1.5:1, 2=2:1 ,3=3:1)
zoom=2              # Zoom (0=none | 1=66% | 2=50% (or 33%)

#---------------------------------------------------------------------
# YOU SHOULD NOT EDIT BELOW
#---------------------------------------------------------------------
# Modify default parameters by input arguments

make=no
case=no
plot_case=no
run=no
data_run=no
plot_run=no
clean=no
all=no

while [ $# -ge 1 ]; do
  case $1 in
    make)        make=yes       ;;
    case)        case=yes       ;;
    plot_case)   plot_case=yes  ;;
    run)         run=yes        ;;
    data_run)    data_run=yes   ;;
    plot_run)    plot_run=yes   ;;
    clean)       clean=yes      ;;
    all)         all=yes        ;;
    check)       check=check    ;;
    outlog)      outlog=outlog  ;;
    -binsrc)     binsrc=$2      ; shift ;;
    -procs)      procs=$2       ; shift ;;
    -vis)        vis=$2         ; shift ;;
    -cordir)     cordir=$2      ; shift ;;
    -heldir)     heldir=$2      ; shift ;;
    -obsdir)     obsdir=$2      ; shift ;;
    -projdir)    projdir=$2     ; shift ;;
    -project)    project=$2     ; shift ;;
    -numo)       numo=$2        ; shift ;;
    -nreg)       nreg=$2        ; shift ;;
    -reg)        reg=$2         ; shift ;;
    -res)        res=$2         ; shift ;;
    -res23)      res23=$2       ; shift ;;
    -npos)       npos=$2        ; shift ;;
    -obsdate)    obsdate=$2     ; shift ;;
    -rundate)    rundate=$2     ; shift ;;
    -tstop)      tstop=$2       ; shift ;;
    -tstep)      tstep=$2       ; shift ;;
    -cordata)    cordata=$2     ; shift ;;
    -corfile)    corfile=$2     ; shift ;;
    -cormode)    cormode=$2     ; shift ;;
    -corobs)     corobs=$2      ; shift ;;
    -corymdh)    corymdh=$2     ; shift ;;
    -amb)        amb=$2         ; shift ;;
    -bndrot)     bndrot=$2      ; shift ;;
    -gamma)      gamma=$2       ; shift ;;
    -xalpha)     xalpha=$2      ; shift ;;
    -cdifb)      cdifb=$2       ; shift ;;
    -qheat)      qheat=$2       ; shift ;;
    -pbnd)       pbnd=$2        ; shift ;;
    -pbnd2)      pbnd2=$2       ; shift ;;
    -pbnd6)      pbnd6=$2       ; shift ;;
    -pbndtim1)   pbndtim1=$2    ; shift ;;
    -pevo1o5)    pevo1o5=$2     ; shift ;;
    -pevo5o1)    pevo5o1=$2     ; shift ;;
    -ptim1)      ptim1=$2       ; shift ;;
    -ptim1o5)    ptim1o5=$2     ; shift ;;
    -ptim2e5)    ptim2e5=$2     ; shift ;;
    -ptim3)      ptim3=$2       ; shift ;;
    -ptim3e1)    ptim3e1=$2     ; shift ;;
    -fontsize)   fontsize=$2    ; shift ;;
    -info)       info=$2        ; shift ;;
    -labitem)    labitem=$2     ; shift ;;
    -rmax)       rmax=$2        ; shift ;;
    -size)       size=$2        ; shift ;;
    -zoom)       zoom=$2        ; shift ;;
  esac
  shift
done

if [ $make == no ] && [ $case == no ] && [ $plot_case == no ] && \
   [ $run == no ] && [ $data_run == no ] && [ $plot_run == no ] && \
   [ $clean == no ] ; then
  make=yes
  case=yes
  plot_case=yes
  run=yes
  data_run=yes
  plot_run=yes
fi

#---------------------------------------------------------------------
# Environment variables

echo source $HELIO_DIR/env/env.inc
     source $HELIO_DIR/env/env.inc

#---------------------------------------------------------------------
# Case codes

casecode=wsa2bc
grdcode=reg2grd

evocode=out2evo
loscode=out2los

#---------------------------------------------------------------------
# Numerical model

# Default values (numo=mp1um1d)
scheme=tvdlf2       # Numerical scheme (tvdlf1 | tvdlf2)
flxdif=noflxdif     # Increased flux diffusion (noflxdif | flxdif)
vsdif=novsdif       # Increased signal speed (novsdif | vsdif)
bchalf=nobchalf     # Time-dependent boundary conditions at the half-step
bctim=bctim1c       # Boundary condition (bctim1c | bctim1h | bctim1v)
difb=difb           # Div(B) diffusion treament (difb | dedner | nodifb)
divb=nodivb         # Div(B) treament (jandif | povell | powell | nodivb)
energy=ethe         # Energy equation (ethe | etot | etot2)
flux=llf            # Numerical flux (llf | hll)
limiter=minmod      # Slope limiter (minmod | albada)
heat=vheat          # Global heating (vheat | hheat | noheat)
mftrace=mfrk4       # IMF line tracing (mfeuler | mfrk4)
upwind=upwind       # Numerical upwinding

# Modification by input argument "numo"
if [ $numo == mp1va2d ] ; then
  upwind=upwindp      # Numerical upwinding
  limiter=albada      # Slope limiter (minmod | albada)
  flux=hll            # Numerical flux (llf | hll)
fi
if [ $numo == mp2um1d ] ; then
  energy=etot         # Energy equation (ethe | etot | etot2)
fi
if [ $numo == mp2va2d ] ; then
  energy=etot         # Energy equation (ethe | etot | etot2)
  upwind=upwindp      # Numerical upwinding
  limiter=albada      # Slope limiter (minmod | albada)
  flux=hll            # Numerical flux (llf | hll)
fi
if [ $numo == mp3um1d ] ; then
  energy=etot2        # Energy equation (ethe | etot | etot2)
fi
if [ $numo == mp3va2d ] ; then
  energy=etot2        # Energy equation (ethe | etot | etot2)
  upwind=upwindp      # Numerical upwinding
  limiter=albada      # Slope limiter (minmod | albada)
  flux=hll            # Numerical flux (llf | hll)
fi

# Label the numerical model by "bchalf"
if [ $bchalf == bchalf ] ; then
  numo=$numo'h'
fi

#---------------------------------------------------------------------
# Computational region and grid parameters

# Default values (reg=1, res=low)
x1l=0.1             # Inner boundary (AU)
x1r=1.1             # Outer boundary (AU)
x2l=30.             # Minimum co-latitude (deg)
x2r=150.            # Maximum co-latitude (deg)
n1=128              # No. of basic X1-grid points
n2=30               # No. of basic X2-grid points
n3=90               # No. of basic X3-grid points

# Modification by input argument "nreg"
if [ $nreg == 2 ] ; then
  n1=144              # No. of basic X1-grid points
fi
if [ $nreg == 3 ] ; then
  n1=96               # No. of basic X1-grid points
fi

# Modification by input arguments "reg" and "nreg"
if [ $nreg != 2 ] ; then
  if [ $reg -gt 1 ] ; then
    x1r=$reg'.1'        # Outer boundary (AU)
    let n1=n1*reg       # No. of basic X1-grid points
  fi
else
  x1r=1.225           # Outer boundary (AU)
  if [ $reg == 2 ] ; then
    x1r=2.35
  fi
  if [ $reg == 3 ] ; then
    x1r=3.475
  fi
  if [ $reg == 4 ] ; then
    x1r=4.6
  fi
  if [ $reg == 5 ] ; then
    x1r=5.725
  fi
  if [ $reg == 6 ] ; then
    x1r=6.85
  fi
  if [ $reg == 7 ] ; then
    x1r=7.975
  fi
  if [ $reg == 8 ] ; then
    x1r=9.1
  fi
  if [ $reg == 9 ] ; then
    x1r=10.225
  fi
  if [ $reg == 10 ] ; then
    x1r=11.35
  fi
  let n1=n1*reg       # No. of basic X1-grid points
fi

# Modification by input argument "res"
nx=1                # Grid resolution factor
if [ $res == med ] ; then
  nx=2                # Grid resolution factor
fi
if [ $res == high ] ; then
  nx=4                # Grid resolution factor
fi
if [ $res == fine ] ; then
  nx=8                # Grid resolution factor
fi
if [ $nx != 1 ] ; then
  let n1=n1*nx
  let n2=n2*nx
  let n3=n3*nx
fi

# Modification by input argument "res23"
if [ $res23 != 1 ] ; then
  let n2=n2*res23     # Refinement of the X2-grid
  let n3=n3*res23     # Refinement of the X3-grid
fi

# Label the grid by "nreg", "reg", "res", "res23", and "npos"
grd=reg
if [ $nreg == 1 ] ; then
  grd=ra
fi
if [ $nreg == 2 ] ; then
  grd=rb
fi
if [ $nreg == 3 ] ; then
  grd=rc
fi
grd=$grd$reg$res'1'
if [ $res23 != 1 ] ; then
  grd=$grd'r'$res23
fi
if [ $npos != 1 ] ; then
  grd=$grd'p'$npos
fi

#---------------------------------------------------------------------
# Computational period and output parameters (relative to "rundate")

# Default values (reg=1)
vutime=3600.        # Time unit scale factor (s)
tstart=-168.        # Start computation at this time (h)
tfrom=0.            # Output from this time (h)

# Modification by input argument "reg"
itstart=${tstart%.*}    # Integer part of "tstart"
let itstart=itstart*reg # Adjust *itstart* for computational domain
tstart=$itstart'.'      # Start computation at this time (h)

# Derived parameters
tfirst=$tstart      # First input data (h)
if [ $cormode == single ] ; then
  tlast=0.          # Last input data (h)
else
  tlast=$tstop      # Last input data (h)
fi
tto=$tstop          # Output to this time (h)
tfstep=$tstep       # FLD output with this step in time (h)
t1step=0.           # SLC1 output with this step in time (h)
t2step=0.           # SLC2 output with this step in time (h)
ttstep=$tstep       # TIM output with this step in time (h)

#---------------------------------------------------------------------
# WSA maps specified by calendar date have syntax by <cordata>:
# (1) <cordata>=wsadu (renamed WSA_v2.* files)renamed by DO and used at CCMC:
# - <cordir>/wsadu/<corobs>/<yyyy>/<mm>/wsadu_<corobs>_<yyyymmddThhmm>.fits
# (2) <cordata>=wsadv (WSA_v2.* files by SWPC):
# - <cordir>/wsadw/<corobs>/<yyyy>/<mm>/<dd>/wsa_vel_21.5rs_<yyyymmddhhmm>_gong.fits
# (3) <cordata>=wsadw (WSA_v4.* files):
# - <cordir>/wsada/<corobs>/<yyyy>/<mm>/vel_<yyyymmddhhmmss>R000_<corobs>.fits
# read by "wsa2bc" code.
#

cordata0=${cordata%/*}  # Get the first directory in cordata

#---------------------------------------------------------------------
# Ambient solar wind parameters

# Default values (amb=a8b1)
bfast=500.          # Radial magnetic field of fast stream (nT)
bslow=400.          # Radial magnetic field of slow stream (nT)
bamb=350.           # Radial magnetic field of mean stream (nT)
if [ ${corobs:0:4} == gong ] ; then
  bscl=5.             # Magnetic field scaling factor
else
  bscl=2.             # Magnetic field scaling factor
fi
dfast=200.          # Number density of fast stream (cm-3)
dslow=4000.         # Number density of slow stream (cm-3)
damb=300.           # Number density of mean stream (cm-3)
dscl=1.             # Number density scaling factor
tfast=1.5           # Mean temperature of fast stream (MK)
tslow=0.1           # Mean temperature of slow stream (MK)
tamb=0.5            # Mean temperature of mean stream (MK)
tscl=1.             # Mean temperature scaling factor
vfast=700.          # Radial flow velocity of fast stream (km/s)
vslow=200.          # Radial flow velocity of slow stream (km/s)
vamb=450.           # Radial flow velocity of mean stream (km/s)
vrfast=25.          # Reduction of the maximum flow velocity (km/s)
vrslow=75.          # Reduction of the minimum flow velocity (km/s)
shift=8.            # Azimuthal shift at the inner boundary (deg)
nshift=1            # Azimuthal shift at the inner boundary (0,1,2,3)
xalpha=0.05         # Fraction of alpha particles (rel. to protons)
dvexp=2.            # Exponent in N*V^dvexp=const condition
nptot=0             # = 0 (1) if P_the (P_tot) balance at boundary
nbrad=1             # Magnetic field correction
nbndace=0           # Adjust WSA values by measured values at ACE (0-3)

# Modification by input argument "amb"
if [ $amb == a3b2 ] ; then
  bfast=300.          # Radial magnetic field of fast stream (nT)
  bslow=0.            # Radial magnetic field of slow stream (nT)
  bamb=0.             # Radial magnetic field of mean stream (nT)
  if [ ${corobs:0:4} == gong ] ; then
    bscl=4.             # Magnetic field scaling factor
  else
    bscl=2.             # Magnetic field scaling factor
  fi
  dfast=200.          # Number density of fast stream (cm-3)
  dslow=2000.         # Number density of slow stream (cm-3)
  damb=0.             # Number density of mean stream (cm-3)
  dscl=1.             # Number density scaling factor
  tfast=0.8           # Mean temperature of fast stream (MK)
  tslow=0.1           # Mean temperature of slow stream (MK)
  tamb=0.             # Mean temperature of mean stream (MK)
  tscl=1.             # Mean temperature scaling factor
  vfast=675.          # Radial flow velocity of fast stream (km/s)
  vslow=225.          # Radial flow velocity of slow stream (km/s)
  vamb=0.             # Radial flow velocity of mean stream (km/s)
  vrfast=25.          # Reduction of the maximum flow velocity (km/s)
  vrslow=25.          # Reduction of the minimum flow velocity (km/s)
  shift=8.            # Azimuthal shift at the inner boundary (deg)
  nshift=1            # Azimuthal shift at the inner boundary (0,1,2,3)
  xalpha=0.05         # Fraction of alpha particles (rel. to protons)
  dvexp=2.            # Exponent in N*V^dvexp=const condition
  nptot=1             # = 0 (1) if P_the (P_tot) balance at boundary
  nbrad=3             # Magnetic field correction
  nbndace=0           # Adjust WSA values by measured values at ACE (0-3)
fi
if [ $amb == a6b1 ] ; then
  bfast=350.          # Radial magnetic field of fast stream (nT)
  bslow=300.          # Radial magnetic field of slow stream (nT)
  bamb=0.             # Radial magnetic field of mean stream (nT)
  if [ ${corobs:0:4} == gong ] ; then
    bscl=4.             # Magnetic field scaling factor
  else
    bscl=2.             # Magnetic field scaling factor
  fi
  dfast=125.          # Number density of fast stream (cm-3)
  dslow=2000.         # Number density of slow stream (cm-3)
  damb=0.             # Number density of mean stream (cm-3)
  dscl=1.             # Number density scaling factor
  tfast=1.5           # Mean temperature of fast stream (MK)
  tslow=0.1           # Mean temperature of slow stream (MK)
  tamb=0.             # Mean temperature of mean stream (MK)
  tscl=1.             # Mean temperature scaling factor
  vfast=700.          # Radial flow velocity of fast stream (km/s)
  vslow=200.          # Radial flow velocity of slow stream (km/s)
  vamb=0.             # Radial flow velocity of mean stream (km/s)
  vrfast=25.          # Reduction of the maximum flow velocity (km/s)
  vrslow=75.          # Reduction of the minimum flow velocity (km/s)
  shift=8.            # Azimuthal shift at the inner boundary (deg)
  nshift=1            # Azimuthal shift at the inner boundary (0,1,2,3)
  xalpha=0.05         # Fraction of alpha particles (rel. to protons)
  dvexp=2.            # Exponent in N*V^dvexp=const condition
  nptot=0             # = 0 (1) if P_the (P_tot) balance at boundary
  nbrad=1             # Magnetic field correction
  nbndace=0           # Adjust WSA values by measured values at ACE (0-3)
fi
if [ $amb == a7b1 ] ; then
  bfast=500.          # Radial magnetic field of fast stream (nT)
  bslow=0.            # Radial magnetic field of slow stream (nT)
  bamb=400.           # Radial magnetic field of mean stream (nT)
  if [ ${corobs:0:4} == gong ] ; then
    bscl=5.             # Magnetic field scaling factor
  else
    bscl=2.             # Magnetic field scaling factor
  fi
  dfast=300.          # Number density of fast stream (cm-3)
  dslow=5000.         # Number density of slow stream (cm-3)
  damb=500.           # Number density of mean stream (cm-3)
  dscl=1.             # Number density scaling factor
  tfast=1.2           # Mean temperature of fast stream (MK)
  tslow=0.1           # Mean temperature of slow stream (MK)
  tamb=0.3            # Mean temperature of mean stream (MK)
  tscl=1.             # Mean temperature scaling factor
  vfast=700.          # Radial flow velocity of fast stream (km/s)
  vslow=200.          # Radial flow velocity of slow stream (km/s)
  vamb=450.           # Radial flow velocity of mean stream (km/s)
  vrfast=20.          # Reduction of the maximum flow velocity (km/s)
  vrslow=50.          # Reduction of the minimum flow velocity (km/s)
  shift=8.            # Azimuthal shift at the inner boundary (deg)
  nshift=1            # Azimuthal shift at the inner boundary (0,1,2,3)
  xalpha=0.05         # Fraction of alpha particles (rel. to protons)
  dvexp=2.            # Exponent in N*V^dvexp=const condition
  nptot=0             # = 0 (1) if P_the (P_tot) balance at boundary
  nbrad=1             # Magnetic field correction
  nbndace=0           # Adjust WSA values by measured values at ACE (0-3)
fi
if [ $amb == ace3a ] ; then
  bfast=500.          # Radial magnetic field of fast stream (nT)
  bslow=0.            # Radial magnetic field of slow stream (nT)
  bamb=0.             # Radial magnetic field of mean stream (nT)
  bscl=0.6            # Magnetic field scaling factor
  dfast=100.          # Number density of fast stream (cm-3)
  dslow=2000.         # Number density of slow stream (cm-3)
  damb=500.           # Number density of mean stream (cm-3)
  dscl=1.             # Number density scaling factor
  tfast=2.            # Mean temperature of fast stream (MK)
  tslow=0.1           # Mean temperature of slow stream (MK)
  tamb=0.6            # Mean temperature of mean stream (MK)
  tscl=1.             # Mean temperature scaling factor
  vfast=700.          # Radial flow velocity of fast stream (km/s)
  vslow=200.          # Radial flow velocity of slow stream (km/s)
  vamb=0.             # Radial flow velocity of mean stream (km/s)
  vrfast=0.           # Reduction of the maximum flow velocity (km/s)
  vrslow=50.          # Reduction of the minimum flow velocity (km/s)
  nbndace=3           # Adjust WSA values by measured values at ACE (0-3)
fi

# Label the ambient solar wind parameters by "gamma"
if [ $gamma == 1.6666667 ] ; then
  amb=$amb'-g53'
else
  zgamma=${gamma//.}
  amb=$amb'-g'${zgamma:0:2}
fi

# Label the ambient solar wind parameters by "xalpha"
if [ $xalpha == 0. ] ; then
  amb=$amb'x00'
else
  zxalpha=${xalpha//.}
  amb=$amb'x'${zxalpha:1:2}
fi

# Label the ambient solar wind parameters by "bndrot"
if [ $bndrot == sidereal ] ; then
  amb=$amb'r'
fi

# Label the coronal data by "cordata"
if [ $cordata0 == wsafr ] ; then
  cor='fr'
else
  if [ $cormode == single ] ; then
    cor='d'
  else
    cor='t'
  fi
  lcorvers=${cordata0: -2}   # Get the last two characters from cordata
  cor=$cor$lcorvers
fi

# Label the coronal data by "corobs"
lcorobs=${corobs: -1}    # Get the last character from corobs
cor=$cor$lcorobs

# Label the boundary conditions
bnd=$cor'-'$amb

#---------------------------------------------------------------------
# Run parameters

# Default values (runpar=g53h10)
ceclip=0.01         # Clip the thermal energy below fraction of the total energy
cflxdif=0.          # Increase flux diffusion (0:1)
cvetot=0.           # Compressive velocity threshold for Etot2 (if Ethe+Etot)
cvsdif=0.           # Increase signal speed (0:1)
nheat=3             # Volumetric heating type
npla=1              # Evolution at planets (0 | 1)
nspa=1              # Evolution at spacecraft (0 | 1)
nfld=3              # Tracing of field lines (0,1,2,3 for none, helio, bnd, helio+bnd)

# Label the run parameters by "nheat" and "qheat"
if [ $nheat == 0 ] ; then
  runpar='q0'
fi
if [ $nheat == 1 ] ; then
  runpar='q'${qheat//.}
fi
if [ $nheat == 2 ] ; then
  runpar='h'${qheat//.}
fi
if [ $nheat == 3 ] ; then
  runpar='t'${qheat//.}
fi

# Label the run parameters by "cvetot"
if [ $energy == etot2 ] ; then
  cvetot=-1.          # Compressive velocity threshold for Etot2 (if Ethe+Etot)
  zzz=${cvetot//.}
  if [ $zzz -lt 0 ] ; then
    zzz='w'${zzz//-}
  else
    zzz='v'$zzz
  fi
  runpar=$runpar$zzz
fi

# Label the run parameters - numerical coefficients
if [ $cdifb != 0. ] ; then
  runpar=$runpar'd'${cdifb//.}
fi

#---------------------------------------------------------------------
# Compose the plot labels

# Plot label - case-observatory
if [ $corobs == gongb ] ; then
 obs=GONGb
else
  if [ $corobs == gongj ] ; then
    obs=GONGj
  else
    if [ $corobs == gongz ] ; then
      obs=GONGz
    else
      obs=${corobs^^}
    fi
  fi
fi

# Plot label - case-data
if [ ${cordata0:0:3} == wsa ] ; then
  wsa='WSA'$lcorvers
  if [ $cormode == single ] ; then
    wsa=$wsa'd'
  else
    wsa=$wsa't'
  fi
else
  wsa=${cordata0^^}
fi

# Plot label - case & run
labcase=$obs'-'$wsa' / '$amb' / '$res
labrun=$obs'-'$wsa' / '$amb' / '$numo' / '$runpar' / '$res

# Plot label - computed at center
labcen='null'
labuser='null'

# Figure name prefix and suffix
namfiga=$bnd
namfigz_case=$res
namfigz_run=$numo'_'$res

#---------------------------------------------------------------------
# Export global parameters

export system batch binsrc vis
export ENLIL_DIR HELIO_DIR
export cordir heldir obsdir
export projdir project numo procs
export casecode grdcode evocode loscode
export grd reg res n1 n2 n3 npos
export obsdate rundate bnd runpar

#---------------------------------------------------------------------
# Compile the codes

if [ $make == yes ] && [ $clean == no ] ; then

  bash $ENLIL_DIR/codes/enlil/scripts/make.sh\
    -check $check -outlog $outlog\
    -energy $energy -heat $heat -bchalf $bchalf -bctim $bctim\
    -scheme $scheme -upwind $upwind -limiter $limiter -flux $flux\
    -flxdif $flxdif -vsdif $vsdif -difb $difb -divb $divb\
    -mftrace $mftrace

fi

#---------------------------------------------------------------------
# Generate boundary conditions

if [ $case == yes ] && [ $clean == no ] ; then

  bash $ENLIL_DIR/codes/enlil/scripts/case.sh\
    -x1l $x1l -x1r $x1r -x2l $x2l -x2r $x2r\
    -vutime $vutime -tfirst $tfirst -tlast $tlast\
    -cordata $cordata -cormode $cormode -corobs $corobs\
    -corymdh $corymdh -corfile $corfile\
    -bndrot $bndrot -shift $shift -nshift $nshift\
    -bfast $bfast -bslow $bslow -bamb $bamb -bscl $bscl\
    -dfast $dfast -dslow $dslow -damb $damb -dscl $dscl\
    -tfast $tfast -tslow $tslow -tamb $tamb -tscl $tscl\
    -vfast $vfast -vslow $vslow -vamb $vamb -vrfast $vrfast -vrslow $vrslow\
    -xalpha $xalpha -dvexp $dvexp -nptot $nptot -nbrad $nbrad\
    -nbndace $nbndace

fi

#---------------------------------------------------------------------
# Plot boundary conditions

if [ $plot_case == yes ] && [ $clean == no ] ; then

  if [ $pbnd != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_case.sh\
      -info $info -labcase "$labcase" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_case"\
      -bcnt 500 -dcnt 4000 -tcnt 2 -vcnt [200,1600]\
      -bevo 500 -devo 2000 -tevo 2 -vevo [200,800]\
      -fontsize $fontsize -size 1 +png\
      -framerate 10 +mp4\
      +plot_bnd
  fi

  if [ $pbnd2 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_case.sh\
      -info $info -labcase "$labcase" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_case"\
      -namobs "earth" -lonplot "obs" +ecliptic\
      -bradcnt 500 -bphicnt 50 -dcnt 4000 -vcnt [200,1600]\
      -heldata $heldata -planets -1 -spacecraft 2\
      +hcs\
      -fontsize $fontsize -size 1 -wide $wide -zoom $zoom +png\
      -framerate 10 +mp4\
      +plot_bnd2
  fi

  if [ $pbnd6 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_case.sh\
      -info $info -labcase "$labcase" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_case"\
      -namobs "earth" -lonplot "geo" +ecliptic\
      -bradcnt 500 -bthecnt 50 -bphicnt 50 -dcnt 4000 -tcnt 2 -vcnt [200,1600]\
      -heldata $heldata -planets -1 -spacecraft 2\
      +hcs\
      -fontsize $fontsize -size 1 -wide $wide -zoom $zoom +png\
      -framerate 10 +mp4\
      +plot_bnd6
  fi

  if [ $pbndtim1 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_case.sh\
      -info $info -labcase "$labcase" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_case"\
      -namcnt ["brad","bphi","btot","den","tem","vrad"]\
      -bcnt 500 -dcnt 4000 -tcnt 2 -vcnt [200,1600]\
      +hcs\
      -fontsize $fontsize -size 1 -wide $wide -zoom $zoom +png\
      +plot_bndtim1
  fi

fi

#---------------------------------------------------------------------
# Run heliospheric computation

if [ $run == yes ] && [ $clean == no ] ; then

  bash $ENLIL_DIR/codes/enlil/scripts/run.sh\
    -vutime $vutime -tstart $tstart -tstop $tstop\
    -tffrom $tfrom -tfto $tto -tfstep $tfstep\
    -t1from $tfrom -t1to $tto -t1step $t1step\
    -t2from $tfrom -t2to $tto -t2step $t2step\
    -ttfrom $tfrom -ttto $tto -ttstep $ttstep\
    -bndrot $bndrot\
    -cdifb $cdifb -cvetot $cvetot -cflxdif $cflxdif -cvsdif $cvsdif\
    -gamma $gamma -xalpha $xalpha\
    -qheat $qheat -nheat $nheat\
    -heldata $heldata -npla $npla -nspa $nspa -nfld $nfld

fi

#---------------------------------------------------------------------
# Process heliospheric computation

if [ $data_run == yes ] && [ $clean == no ] ; then

  bash $ENLIL_DIR/vis/idl-gdl/scripts/data_run.sh\
    -vutime $vutime -xtmin 0. -xtmax $tstop\
    +data_evo

fi

#---------------------------------------------------------------------
# Plot heliospheric computation

if [ $plot_run == yes ] && [ $clean == no ] ; then

  if [ $pevo1o5 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth","stereoa"]\
      -namevo ["vr"] -vrevo [200,800]\
      +evohcld +evopos +evosho -xpsho 1.2\
      -vutime $vutime -xtobs [$tfrom,$tto] -xtplot [$tfrom,$tto]\
      -fontsize $fontsize -size 1 -wide $wide +png\
      +plot_evo1o5
  fi

  if [ $pevo5o1 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth"]\
      -namevo ["vr","den","tem","btot"] -vrevo [200,800] -devo 40 -tevo 1.0 -bevo 20\
      +evopos\
      -vutime $vutime -xtobs [$tfrom,$tto] -xtplot [$tfrom,$tto]\
      -fontsize $fontsize -size 1 -wide $wide +png\
      +plot_evo5o1
  fi

  if [ $ptim1 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth"] -lonplot "geo"\
      -namfld ["earth","stereoa"]\
      -namcnt ["vrad","den"] -vcnt [200,1600] -dcnt 60\
      +hcs -rmax $rmax -planets -1 -spacecraft 2 +ecliptic\
      -vutime $vutime -xtplot [$tfrom,$tto] -istep $ptim1\
      -fontsize $fontsize -size 1 -wide $wide -zoom $zoom +png\
      -framerate 20 +mp4\
      +plot_tim1
  fi

  if [ $ptim1o5 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth","stereoa"] -lonplot "geo"\
      -namfld ["earth","stereoa"]\
      -namcnt ["vrad"] -vcnt [200,1600] -dcnt 60\
      -namevo ["vr"] -vrevo [200,800]\
      +evopos\
      +hcs -rmax $rmax -planets -1 -spacecraft 2 +ecliptic\
      -vutime $vutime -xtobs [$tfrom,$tto] -xtplot [$tfrom,$tto] -istep $ptim1o5\
      -fontsize $fontsize -size 1 -zoom $zoom +png\
      -framerate 20 +mp4\
      +plot_tim1o5
  fi

  if [ $ptim2e5 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth"] -lonplot "geo"\
      -namfld ["earth","stereoa"]\
      -namcnt ["vrad","den"] -vcnt [200,1600] -dcnt 60\
      -namevo ["vr","den","tem","btot"] -vrevo [200,800] -devo 60\
      +evopos\
      +hcs -rmax $rmax -planets -1 -spacecraft 2 +ecliptic\
      -vutime $vutime -xtobs [$tfrom,$tto] -xtplot [$tfrom,$tto] -istep $ptim2e5\
      -fontsize $fontsize -size 1 -zoom $zoom +png\
      -framerate 20 +mp4\
      +plot_tim2e5
  fi

  if [ $ptim3e1 != 0 ] ; then
    bash $ENLIL_DIR/vis/idl-gdl/scripts/plot_run.sh\
      -info $info -labrun "$labrun" -labcen "$labcen" -labuser "$labuser"\
      -labitem "$labitem" -namfiga "$namfiga" -namfigz "$namfigz_run"\
      -namobs ["earth"] -lonplot "geo"\
      -namfld ["earth","stereoa"]\
      -namcnt ["vrad","den"] -vcnt [200,1600] -dcnt 60\
      -vrevo [200,800] -devo 60\
      +evopos\
      +hcs -rmax $rmax -planets -1 -spacecraft 2 +ecliptic\
      -vutime $vutime -xtobs [$tfrom,$tto] -xtplot [$tfrom,$tto] -istep $ptim3e1\
      -fontsize $fontsize -size 1 -zoom $zoom +png\
      -framerate 20 +mp4\
      +plot_tim3e1
  fi

fi

#---------------------------------------------------------------------
# Remove non-critical files

if [ $clean == yes ] ; then

  bash $ENLIL_DIR/codes/enlil/scripts/clean.sh\
    -make $make -case $case -run $run -all $all

fi

#---------------------------------------------------------------------
exit 0