#!/bin/bash
#=====================================================================
#
# NAME 
#   xwambfr
#
# PURPOSE
#   Master script to run *ENLIL* with *WSA* full-rotation 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> = wsafr<vers> indicates
#   data from different WSA versions and <corobs> = gong, nso
# 
#   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
#   xwambfr [-<parameter> <value>]
#   xwambfr make [-<parameter> <value>]
#   xwambfr case [-<parameter> <value>]
#   xwambfr plot_case [-<parameter> <value>]
#   xwambfr run [-<parameter> <value>]
#   xwambfr data_run [-<parameter> <value>]
#   xwambfr plot_run [-<parameter> <value>]
#   xwambfr 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/wambfr
#
# 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:
#                       2010-03-26 | <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:
#                       wsafr
#   -corfile <char>   - Coronal data file (null if cordata used):
#                       null | <name>
#   -cormode <char>   - Coronal data mode:
#                       single | multi
#   -corobs <char>    - Solar photospheric field observatory:
#                       nso | gongb | gongz | mdi | mwo
#   -crnum <char>     - Carrington Rotation number:
#                       2095 | <xxxx>
#   -crlon <char>     - Carrington Rotation leading longitude:
#                       null | <xxx>
#   -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/wambfr # 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=2010-03-26  # Start computation at this date (null | yyyy-mm-ddThh)
tstop=672.          # Stop computation at this time (h)
tstep=336.          # Output with this step in time (h)

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

# Input data
cordata=wsafr       # Coronal data sub-directory (wsafr)
corfile=null        # Coronal data file (null if cordata used)
cormode=single      # Coronal data mode (single | multi)
corobs=nso          # Observatory name (gongb | gongz | mdi | mwo | nso)
crnum=2095          # Carrington Rotation number (null | xxxx)
crlon=null          # Carrington Rotation leading longitude (null | xxx)

# 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

info=1              # Informative footnote (1=yes, 0=no, -1=ghost)
fontsize=1.0        # Font size (1=small, 2=medium, 3=large)
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 ;;
    -crnum)      crnum=$2       ; shift ;;
    -crlon)      crlon=$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)

#---------------------------------------------------------------------
# Coronal data sub-directory for a WSA maps specified by Carrington
# rotation with a syntax (original WSA approach):
# - <cordir>/<cordata>/vel_<crnum>_<crlon>.00_01_<corobs>.fits
# read by "wsa2bc" code.
# NOTE: WSA full-rotation maps have <crlon>=000.
#
# WSA maps specified by calendar date have syntax by <cordata>:
# (1) <cordata>=wsafr (renamed WSA_v2.* files)renamed by DO and used at CCMC:
# - <cordir>/wsafr/<corobs>/<yyyy>/<mm>/wsadu_<corobs>_<yyyymmddThhmm>.fits
# read by "wsa2bc" code.
#

cordata0=${cordata%/*}  # Get the first directory in cordata
if [ $corobs == gong ] ; then
  cordata=$cordata'/WSA_VEL/GONG/CARR/2.5deg/L72R21.5'
fi
if [ $corobs == gongb ] ; then
  cordata=$cordata'/WSA_VEL/GONG/CARR/2.5deg/L72R21.5'
fi
if [ $corobs == mdi ] ; then
  cordata=$cordata'/WSA_VEL/MDI/CARR/2.5deg/L72R21.5'
fi
if [ $corobs == mwo ] ; then
  cordata=$cordata'/WSA_VEL/MWO/CARR/2.5deg/L72R21.5'
fi
if [ $corobs == nso ] ; then
  cordata=$cordata'/WSA_VEL/NSO/CARR/2.5deg/L72R21.5'
fi
#if [ $cordata0 == wsafr ] ; then
  corfile='vel_'$crnum'_000.00_01_'$corobs'.fits'
#fi

#---------------------------------------------------------------------
# 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=WSAfr
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 crnum crlon
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\
    -cordir $PWD -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\
      +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 1 +png\
      +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 1 +png\
      +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 1 +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 20 -tevo 60 -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