Commits

Tom Roche  committed 9286789

make hourly output

Added make_hourlies.ncl and integrated into drivers (AQMEII_driver.sh and uber_driver.sh) to
* create netCDF output with |TSTEP|=25
* move NCL-compliant netCDF to *.ncf
* launch VERDI on *.ncf

Also tweaks to distill_EPIC_emittivities.ncl (from which make_hourlies.ncl is derived).

Tested on terrae/EMVL from both AQMEII_driver.sh and uber_driver.sh::file (latter in fresh terminal), and VERDI opens output.

TODO:
* handle {name, path to} VERDI in regrid_utils/bash_utilities.sh
* generate plot of EPIC-EDGAR (i.e., show where EDGAR < EPIC, EDGAR > EPIC)
*** add to README, wiki: links to overviews of BELD, EPIC
*** add to wiki: Simulation-of-N2O-Production-and-Transport-in-the-US-Cornbelt-Compared-to-Tower-Measurements#wiki-input-processing-EPIC-0509 (I have input processing instructions for other EIs, but not EPIC)
* add more stats to README, esp comparing before EPIC overlay to after
* move explicit constants BELD_CSV_to_RDS.r -> AQMEII_driver.sh
* regrid_utils: make commandline netCDF.stats.to.stdout.r take files beginning with numerals, e.g.
- $ Rscript ./netCDF.stats.to.stdout.r netcdf.fp=5yravg_20111219_pure_emis.nc data.var.name=N2O # fails
+ $ Rscript ./netCDF.stats.to.stdout.r netcdf.fp=yravg_20111219_pure_emis.nc data.var.name=N2O # works
* port AQMEII_driver failfast to
** other *_driver: CLMCN_driver, EDGAR_driver, GEIA_driver, GFED_driver
** other uber_driver: GFED::uber_driver
* create common project for regrid_resources à la regrid_utils, so I don't hafta hunt down which resource is in which project
* (this project only) use the EPIC input as template, rather than the TEMPLATE_INPUT_FP I shlep to all the projects.
* all regrids: how to nudge off/onshore as required? e.g., soil or burning emissions should never be offshore, marine emissions should never be onshore.
* all regrid maps: add Caribbean islands (esp Bahamas! for offshore burning), Canadian provinces, Mexican states
* complain to ncl-talk about NCL "unsupported extensions," e.g., .ncf and <null/> (e.g., MCIP output)
* determine why '<-' assignment is occasionally required in calls to visualize.*(...)
* document platform versions (e.g., linux, NCL, R)
* test on
** tlrPanP5 (which now has R package=ncdf4, but readAsciiTable of input .txt's is very slow compared to terrae)
** HPCC (once problem with ncdf4 on amad1 is debugged: in process with JOB and KMF)

  • Participants
  • Parent commits 253ab91

Comments (0)

Files changed (4)

File AQMEII_driver.sh

 export COMPUTE_EPIC_FP="${WORK_DIR}/${COMPUTE_EPIC_FN}"
 export COMBINE_EMIS_FN='combine_EDGAR_and_EPIC_emissions.r'
 export COMBINE_EMIS_FP="${WORK_DIR}/${COMBINE_EMIS_FN}"
-# ----------------------------------------------------------------------
 export MAKE_HOURLIES_FN='make_hourlies.ncl'
 export MAKE_HOURLIES_FP="${WORK_DIR}/${MAKE_HOURLIES_FN}"
 export CHECK_CONSERVATION_FN='check_conservation.ncl'
 export AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST="${AQMEIINA_DV_LONG_NAME}"
 export AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST_TYPE='text' # or 'character': ncdf4 types
 # IOAPI pads fileattr=filedesc to length=60 with trailing spaces
-export AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC="$(printf '%-16s' 'EPIC+EDGAR N2O mol/s/yr')"
+export AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC="$(printf '%-16s' 'EPIC+EDGAR N2O mol/s')"
 export AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC_TYPE='character'
 
 ## visualization
 export AQMEII_BOTH_N2O_EMISS_PDF_FP="${WORK_DIR}/${AQMEII_BOTH_N2O_EMISS_PDF_FN}"
 
 # ----------------------------------------------------------------------
-# final output: hourly emissions over AQMEII
+# final output: hourly combined emissions over AQMEII
 # ----------------------------------------------------------------------
 
-### template for multiple netCDF output files (one per day)
-export GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING='@@@@'
-# CMAQ wants an extension that NCL won't write
-GFED_REGRID_HOURLY_EMIS_FN_TEMPLATE_CMAQ="emis_mole_N2O_${MODEL_YEAR}${GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING}.${NETCDF_EXT_CMAQ}"
-GFED_REGRID_HOURLY_EMIS_FN_TEMPLATE_NCL="${GFED_REGRID_HOURLY_EMIS_FN_TEMPLATE_CMAQ}.${NETCDF_EXT_NCL}"
-export GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_CMAQ="${WORK_DIR}/${GFED_REGRID_HOURLY_EMIS_FN_TEMPLATE_CMAQ}"
-export GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL="${WORK_DIR}/${GFED_REGRID_HOURLY_EMIS_FN_TEMPLATE_NCL}"
+# We output annual values, so make output for no specific date
+AQMEII_HOURLY_N2O_EMISS_ROOT="emis_mole_N2O_${MODEL_YEAR}_12US1_cmaq_cb05_soa_2008ab_08c"
+AQMEII_HOURLY_N2O_EMISS_FN_NCL="${AQMEII_HOURLY_N2O_EMISS_ROOT}.${NETCDF_EXT_NCL}"
+export AQMEII_HOURLY_N2O_EMISS_FN_CMAQ="${AQMEII_HOURLY_N2O_EMISS_ROOT}.${NETCDF_EXT_CMAQ}"
+# generally I'll work with the NCL-friendly filename ... because I hafta :-(
+export AQMEII_HOURLY_N2O_EMISS_FP_NCL="${WORK_DIR}/${AQMEII_HOURLY_N2O_EMISS_FN_NCL}"
+export AQMEII_HOURLY_N2O_EMISS_FP="${AQMEII_HOURLY_N2O_EMISS_FN_NCL}"
+# ... but this driver will create this ultimate output
+export AQMEII_HOURLY_N2O_EMISS_FP_CMAQ="${WORK_DIR}/${AQMEII_HOURLY_N2O_EMISS_FN_CMAQ}"
 
-## just the 3-hour regrids for January 01
-export GFED_REGRID_HOURLY_EMIS_JAN_01_FP_CMAQ="${GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_CMAQ}/${GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING}/0101}"
-export GFED_REGRID_HOURLY_EMIS_JAN_01_FP_NCL="${GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL}/${GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING}/0101}"
+## dimensions
+export AQMEII_HOURLY_N2O_EMISS_DIM_TSTEP_N='25'
+# names: unfortunately ncdf4 appears to need these defined to walk its ADSs :-(
+export AQMEII_HOURLY_N2O_EMISS_DIM_X_NAME="${AQMEIINA_DIM_X_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_DIM_Y_NAME="${AQMEIINA_DIM_Y_NAME}"
+# might as well make these explicit, if we're gonna create coordvars
+export AQMEII_HOURLY_N2O_EMISS_DIM_X_ATTR_UNITS="${AQMEIINA_DIM_X_UNITS}"
+export AQMEII_HOURLY_N2O_EMISS_DIM_X_ATTR_LONG_NAME="${AQMEIINA_DIM_X_LONG_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_DIM_Y_ATTR_UNITS="${AQMEIINA_DIM_Y_UNITS}"
+export AQMEII_HOURLY_N2O_EMISS_DIM_Y_ATTR_LONG_NAME="${AQMEIINA_DIM_Y_LONG_NAME}"
 
-## file/global attributes
-export GFED_REGRID_HOURLY_EMIS_HISTORY="see ${PROJECT_WEBSITE}"
-export GFED_REGRID_HOURLY_EMIS_TITLE="hourly N2O emissions due to biomass burning for ${MODEL_YEAR}, derived from GFED emissions and fractions"
+## datavar
+export AQMEII_HOURLY_N2O_EMISS_DATAVAR_NAME="${AQMEIINA_DV_NAME}" # dim=(ROW, COL)
+export AQMEII_HOURLY_N2O_EMISS_DATAVAR_ATTR_LONG_NAME="${AQMEIINA_DV_LONG_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_DATAVAR_ATTR_UNITS="${AQMEIINA_DV_UNITS}"
+export AQMEII_HOURLY_N2O_EMISS_DATAVAR_ATTR_VAR_DESC="${AQMEIINA_DV_VAR_DESC}"
+# "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
+export AQMEII_HOURLY_N2O_EMISS_TFLAG_NAME="${AQMEIINA_TFLAG_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_TFLAG_ATTR_LONG_NAME="${AQMEIINA_TFLAG_LONG_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_TFLAG_ATTR_UNITS="${AQMEIINA_TFLAG_UNITS}"
+export AQMEII_HOURLY_N2O_EMISS_TFLAG_ATTR_VAR_DESC="${AQMEIINA_TFLAG_VAR_DESC}"
 
-### datavar
-export GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME='N2O'
-## format datavar attributes IOAPI-style
-# IOAPI pads varattr=units to length=16 with trailing spaces
-export GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_UNITS="$(printf '%-16s' 'moles/s')"
-# IOAPI pads varattr=long_name to length=16 with trailing spaces
-export GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_LONG_NAME="$(printf '%-16s' ${GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME})"
-# IOAPI pads varattr=var_desc to length=80 with trailing spaces
-export GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_VAR_DESC="$(printf '%-80s' ${GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME})"
+## file/global attributes
+# copy these from input
+export AQMEII_HOURLY_N2O_EMISS_ATTR_NLAYS='1' # output semantics: "layer" == vertical, not crop
+export AQMEII_HOURLY_N2O_EMISS_ATTR_NLAYS_TYPE='integer' # 'short' -> `:NLAYS = 1s ;`
+export AQMEII_HOURLY_N2O_EMISS_ATTR_VARLIST="${AQMEIINA_DV_LONG_NAME}"
+export AQMEII_HOURLY_N2O_EMISS_ATTR_VARLIST_TYPE='text' # or 'character': ncdf4 types
+# IOAPI pads fileattr=filedesc to length=60 with trailing spaces
+export AQMEII_HOURLY_N2O_EMISS_ATTR_FILEDESC="${AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC}"
+export AQMEII_HOURLY_N2O_EMISS_ATTR_FILEDESC_TYPE="${AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC_TYPE}"
 
-# datavar dims=(TSTEP, LAY, ROW, COL)
-export GFED_REGRID_HOURLY_EMIS_N_TSTEP="$(( HOURS_PER_DAY+1 ))"
-export GFED_REGRID_HOURLY_EMIS_N_LAYER='1' # one species in file
-
-# # datavar dims and attributes
-# export GFED_REGRID_HOURLY_EMIS_DIM_TSTEP_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TSTEP_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_TSTEP_UNITS="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TSTEP_UNITS}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_TSTEP_LONG_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TSTEP_LONG_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_LAYER_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_LAYER_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_LAYER_LONG_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_LAYER_LONG_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_LAYER_UNITS="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_LAYER_UNITS}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_X_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_X_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_X_UNITS="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_X_UNITS}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_X_LONG_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_X_LONG_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_Y_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_Y_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_Y_UNITS="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_Y_UNITS}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_Y_LONG_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_Y_LONG_NAME}"
+## no visualization: will do some in check_conservation
 
 # ----------------------------------------------------------------------
 # functions
 #  echo -e "\n$ ${THIS_FN}: PDF_VIEWER='${PDF_VIEWER}'" # debugging
 } # end function setup
 
-#     'get_check_conservation' \
-#     'get_process_EPIC' \
-#     'get_overlay_EPIC' \
-#     'get_make_hourlies' \
+#    'get_check_conservation' \
 function get_helpers {
   for CMD in \
     'get_regrid_utils' \
     'get_string_funcs' \
     'get_compute_EPIC' \
     'get_combine_emis' \
+    'get_make_hourlies' \
   ; do
     echo -e "\n$ ${THIS_FN}::${FUNCNAME[0]}::${CMD}\n"
     eval "${CMD}" # comment this out for NOPing, e.g., to `source`
   fi
 } # end function get_combine_emis
 
+function get_make_hourlies {
+  if [[ -z "${MAKE_HOURLIES_FP}" ]] ; then
+    echo -e "${THIS_FN}: ERROR: MAKE_HOURLIES_FP not defined"
+    exit 1
+  fi
+  # is in this repo
+#  if [[ ! -r "${MAKE_HOURLIES_FP}" ]] ; then
+#    for CMD in \
+#      "${WGET_TO_FILE} ${MAKE_HOURLIES_FP} ${MAKE_HOURLIES_URI}" \
+#    ; do
+#      echo -e "$ ${CMD}"
+#      eval "${CMD}"
+#    done
+#  fi
+  if [[ ! -r "${MAKE_HOURLIES_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: MAKE_HOURLIES_FP=='${MAKE_HOURLIES_FP}' not readable"
+    exit 1
+  fi
+} # end function get_make_hourlies
+
+function get_check_conservation {
+  if [[ -z "${CHECK_CONSERVATION_FP}" ]] ; then
+    echo -e "${THIS_FN}: ERROR: CHECK_CONSERVATION_FP not defined"
+    exit 1
+  fi
+  # is in this repo
+#  if [[ ! -r "${CHECK_CONSERVATION_FP}" ]] ; then
+#    for CMD in \
+#      "${WGET_TO_FILE} ${CHECK_CONSERVATION_FP} ${CHECK_CONSERVATION_URI}" \
+#    ; do
+#      echo -e "$ ${CMD}"
+#      eval "${CMD}"
+#    done
+#  fi
+  if [[ ! -r "${CHECK_CONSERVATION_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: CHECK_CONSERVATION_FP=='${CHECK_CONSERVATION_FP}' not readable"
+    exit 1
+  fi
+} # end function get_check_conservation
+
 function setup_resources {
   for CMD in \
     'get_EDGAR_N2O' \
   done
 } # end function BELD_CSV_to_RDS
 
-### From raw EPIC, "distill" or purify only the N2O emittivities, and demononticize them if necessary.
+### From raw EPIC, "distill" or purify only the N2O emittivities, and demonotonicize them if necessary, and reunit them.
 ### Depends on function get_distill_EPIC above.
 ### Assumes `ncl` is in PATH: see bash_utilities.sh
 ### TODO: pass commandline args to NCL
 
 } # end function combine_emis
 
+### Create (even more) CMAQ-friendly 25-hour *.ncf files.
+### Depends on function get_make_hourlies above.
+### Assumes `ncl` is in PATH: see bash_utilities.sh
+### TODO: pass commandline args to NCL
+### punt: just use envvars :-(
+function make_hourlies {
+  if [[ -z "${AQMEII_HOURLY_N2O_EMISS_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: AQMEII_HOURLY_N2O_EMISS_FP not defined"
+    exit 123
+  fi
+  if [[ -z "${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: AQMEII_HOURLY_N2O_EMISS_FP_CMAQ not defined"
+    exit 124
+  fi
+  if [[ -r "${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: skipping ${MAKE_HOURLIES_FN}: output='${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}' is readable"
+  else # make it ...
+
+    # ... for which we need our input
+    if [[ -z "${AQMEII_BOTH_N2O_EMISS_FP}" ]] ; then
+      echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: AQMEII_BOTH_N2O_EMISS_FP not defined"
+      exit 125
+    fi
+
+    if [[ ! -r "${AQMEII_BOTH_N2O_EMISS_FP}" ]] ; then
+      echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: raw EPIC emittivities=='${AQMEII_BOTH_N2O_EMISS_FP}' not readable"
+      exit 126
+    fi
+
+    if [[ -z "${MAKE_HOURLIES_FP}" ]] ; then
+      echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: MAKE_HOURLIES_FP not defined"
+      exit 127
+    fi
+
+    if [[ ! -r "${MAKE_HOURLIES_FP}" ]] ; then
+      echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: MAKE_HOURLIES_FP=='${MAKE_HOURLIES_FP}' not readable"
+      exit 128
+    fi
+
+#     ncl # bail to NCL and copy script lines, or ...
+    for CMD in \
+      "ncl -n ${MAKE_HOURLIES_FP}" \
+      "mv ${AQMEII_HOURLY_N2O_EMISS_FP} ${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" \
+    ; do
+      cat <<EOM
+
+About to run command="${CMD}"
+
+EOM
+      eval "${CMD}"
+      if [[ $? -ne 0 ]] ; then
+        echo -e "$ ${THIS_FN}::${FUNCNAME[0]}::${CMD}: ERROR: failed or not found\n"
+        exit 129
+      fi
+    done
+  fi # end if [[ -r "${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" ]] ; then
+
+  if [[ -r "${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" ]] ; then
+    # get/run VERDI command here rather than NCL, because we move the NCL-happy output before the NCL-launched process can get to it
+    AQMEII_HOURLY_N2O_EMISS_FP_CMAQ_FQN="$(readlink -f ${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ})" # VERDI really, really wants absolute paths
+    AQMEII_HOURLY_N2O_EMISS_DATAVAR_NAME_VERDI="${AQMEII_HOURLY_N2O_EMISS_DATAVAR_NAME}[1]"
+    AQMEII_HOURLY_N2O_EMISS_VERDI_CMD="verdi -f ${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ_FQN} -s ${AQMEII_HOURLY_N2O_EMISS_DATAVAR_NAME_VERDI} -gtype tile &"
+
+    for CMD in \
+      "ls -alt ${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" \
+      "ncdump -h ${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}" \
+      "${AQMEII_HOURLY_N2O_EMISS_VERDI_CMD}" \
+    ; do
+      echo -e "\n$ ${THIS_FN}::${FUNCNAME[0]}::${CMD}\n"
+      eval "${CMD}" # comment this out for NOPing, e.g., to `source`
+      if [[ $? -ne 0 ]] ; then
+        echo -e "$ ${THIS_FN}::${FUNCNAME[0]}::${CMD}: ERROR: failed or not found\n"
+        exit 130
+      fi
+    done
+  else
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: could not create hourly combined emissions=='${AQMEII_HOURLY_N2O_EMISS_FP_CMAQ}'"
+    exit 131
+  fi
+
+} # end function make_hourlies
+
 function teardown {
 #  echo -e 'TODO!'
   for CMD in \
 #   'distill_EPIC' \          # Separate EPIC N2O emittivities and demonotonicize
 #   'compute_EPIC' \          # Compute EPIC N2O emissions, write, and visualize
 #   'combine_emis' \          # EPIC values replace EDGAR where both exist
-# ----------------------------------------------------------------------
-#   'check_conservation' \    # of mass from inputs to outputs, and check EPIC vs EDGAR values
 #   'make_hourlies' \         # create our real output: hourly emissions for CMAQ
+#   'check_conservation' \    # of mass from inputs to outputs, and check EPIC vs EDGAR values
 #   'teardown' \
 # should always
 # * begin with `setup` to do `module add`
   'distill_EPIC' \
   'compute_EPIC' \
   'combine_emis' \
+  'make_hourlies' \
   'teardown' \
 ; do
   echo -e "\n$ ${THIS_FN}::main loop::${CMD}\n"

File distill_EPIC_emittivities.ncl

 ;!/usr/bin/env ncl ; requires version >= 6.0.0 for default_fillvalue(...)
 ;----------------------------------------------------------------------
-;;; TODO: describe me!
+;;; "distill" raw EPIC emittivities by
+;;; * removing non-N2O emittivities
+;;; * demonotonicize them (if necessary)
+;;; * reunit-ing them
 
 ;----------------------------------------------------------------------
 ; libraries
   pure_emitt_dv_name_VERDI = pure_emitt_dv_name+"[1]"
   ;; VERDI really wants an FQP
   pure_emitt_fp_VERDI = systemfunc("readlink -f "+pure_emitt_fp)
-  pure_emitt_dv_VERDI_cmd = "verdi -f "+pure_emitt_fp_VERDI+" -s "+pure_emitt_dv_name_VERDI+" -gtype tile"
+  pure_emitt_dv_VERDI_cmd = "verdi -f "+pure_emitt_fp_VERDI+" -s "+pure_emitt_dv_name_VERDI+" -gtype tile &"
 
 ;   print(str_get_nl() + \
 ;     this_fn+"; debugging IOAPI: about to run '"+pure_emitt_dv_VERDI_cmd+"'" + str_get_nl())

File make_hourlies.ncl

+;!/usr/bin/env ncl ; requires version >= ???
+;----------------------------------------------------------------------
+; Create CMAQ-friendly 25-hour *.ncf files. Same units (mol/s) as previous.
+
+;----------------------------------------------------------------------
+; libraries
+;----------------------------------------------------------------------
+
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"    ; all built-ins?
+;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; for int2flt,  *_Var*
+
+;----------------------------------------------------------------------
+; constants
+;----------------------------------------------------------------------
+
+this_fn = getenv("MAKE_HOURLIES_FN") ; for debugging
+
+;;; combined (and regridded and reunited) emissions are our input
+input_fp = getenv("AQMEII_BOTH_N2O_EMISS_FP")                  ; path
+input_dv_name = getenv("AQMEII_BOTH_N2O_EMISS_DATAVAR_NAME")   ; real data
+input_TFLAG_name = getenv("AQMEII_BOTH_N2O_EMISS_TFLAG_NAME")  ; IOAPI TFLAG
+
+;;; CMAQ-friendly 25-hour output, but not *.ncf (driver handles that, since NCL won't)
+output_fp = getenv("AQMEII_HOURLY_N2O_EMISS_FP_NCL")
+output_dv_name = getenv("AQMEII_HOURLY_N2O_EMISS_DATAVAR_NAME")
+output_TFLAG_name = getenv("AQMEII_HOURLY_N2O_EMISS_TFLAG_NAME")
+output_TFLAG_tsteps_n = stringtoint(getenv("AQMEII_HOURLY_N2O_EMISS_DIM_TSTEP_N"))
+output_TFLAG_day = 2008001
+
+;----------------------------------------------------------------------
+; functions
+;----------------------------------------------------------------------
+
+;----------------------------------------------------------------------
+; code
+;----------------------------------------------------------------------
+
+begin ; skip if copy/paste-ing to console
+
+;----------------------------------------------------------------------
+; payload
+;----------------------------------------------------------------------
+
+;----------------------------------------------------------------------
+; load input datavar
+;----------------------------------------------------------------------
+
+  input_fh = addfile(input_fp, "r")
+  input_dv = input_fh->$input_dv_name$
+  input_dv_dims_sizes = dimsizes(input_dv)
+  input_dv_dims_n = dimsizes(input_dv_dims_sizes) ; == |dims|
+  input_dv_dims_names = getvardims(input_dv)
+  input_dv_tsteps_n = input_dv_dims_sizes(0)
+  input_dv_layers_n = input_dv_dims_sizes(1) ; LAY == vertical
+  input_dv_rows_n = input_dv_dims_sizes(2)
+  input_dv_cols_n = input_dv_dims_sizes(3)
+  input_dv_mv = getFillValue(input_dv) ; missing value
+  input_dv_type = typeof(input_dv)     ; type as string
+;   input_dv_units = input_dv@units
+
+;----------------------------------------------------------------------
+; load input TFLAG
+;----------------------------------------------------------------------
+
+  ;;; For IOAPI compatibility (mostly so EJC can load this into VERDI), we need IOAPI's fake "datavar" name=TFLAG
+  input_TFLAG = input_fh->TFLAG
+  input_TFLAG_type = typeof(input_TFLAG)     ; type as string
+  input_TFLAG_mv = getFillValue(input_TFLAG) ; missing value
+
+;;; Major aside: IOAPI does not apparently typically set attr={_FillValue, missing_value} on datavar=TFLAG, which makes input_TFLAG_mv==NA
+;;; Below we then propagate that to pure_emitt_TFLAG, and NCL is good with that.
+;;; But when we try to use this with R::ncdf4, it causes problems. So set input_TFLAG_mv, and
+;;; TODO: raise issue on netcdfgroup@unidata.ucar.edu: NCL OK with putting datavar with missing _FillValue, ncdf4 is not
+
+  if (ismissing(input_TFLAG_mv)) then
+    input_TFLAG_mv = default_fillvalue(input_TFLAG_type)
+  end if
+  input_TFLAG_dims_sizes = dimsizes(input_TFLAG)
+  input_TFLAG_dims_n = dimsizes(input_TFLAG_dims_sizes) ; == |dims|
+  input_TFLAG_dims_names = getvardims(input_TFLAG) ; TSTEP, VAR, DATE-TIME
+
+;   ; start debugging----------------------------------------------------
+;   print("printVarSummary(input_TFLAG)==")
+;   printVarSummary(input_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+;----------------------------------------------------------------------
+; setup output datavar with more TSTEP
+;----------------------------------------------------------------------
+
+  ;;; TODO: should be doing this @ driver level
+  output_dv_cols_n = input_dv_cols_n
+  output_dv_layers_n = input_dv_layers_n
+  output_dv_rows_n = input_dv_rows_n
+  output_dv_tsteps_n = output_TFLAG_tsteps_n
+  output_dv_type = input_dv_type
+  output_dv_mv = input_dv_mv
+  output_dv_dims_names = input_dv_dims_names
+  output_dv_dims_sizes = \
+    (/ output_dv_tsteps_n, output_dv_layers_n, output_dv_rows_n, output_dv_cols_n /)
+
+  ;;; output data is same as input data (mol/s over grid), just more of it (more TSTEPs)
+  output_arr = new(output_dv_dims_sizes, output_dv_type, output_dv_mv)
+
+  ;;; iterate dims=(TSTEP, LAY, ROW, COL)
+;   do i_tstep = 0, 0 ; is all we're gonna do, but ...
+  do i_tstep = 0, output_dv_tsteps_n-1
+;     do i_layer = 0, 0
+    do i_layer = 0, output_dv_layers_n-1
+      ; copy grid from the (ASSERT) only layer and tstep in the input
+      output_arr(i_tstep, i_layer,:,:) = (/ input_dv(0,0,:,:) /)
+    end do
+  end do
+
+;----------------------------------------------------------------------
+; setup output TFLAG with more TSTEP
+;----------------------------------------------------------------------
+
+  ;;; output should have
+  ;;; * more |TSTEP|: output_TFLAG_tsteps_n set above
+  ;;; * same |VAR| = 1 (only one "real" datavar=N2O)
+  output_TFLAG_vars_n = input_TFLAG_dims_sizes(1)
+  ;;; * same |DATE-TIME| (constant == 2 for IOAPI)
+  output_TFLAG_datetimes_n = input_TFLAG_dims_sizes(2)
+  output_TFLAG_dims_sizes = \
+    (/ output_TFLAG_tsteps_n, output_TFLAG_vars_n, output_TFLAG_datetimes_n /)
+  output_TFLAG_type = input_TFLAG_type
+  output_TFLAG_dims_names = input_TFLAG_dims_names
+  output_TFLAG_mv = input_TFLAG_mv
+
+  ;;; make the "datavar", and ...
+  output_TFLAG = \
+    new(output_TFLAG_dims_sizes, output_TFLAG_type, output_TFLAG_mv)
+  ;;; ... copy its dimensions and attributes
+;  copy_VarMeta(input_TFLAG, output_TFLAG) ; our dimensions are not a "leftmost subset"
+  do i_dim = 0 , input_TFLAG_dims_n-1
+    ; note subscripting/indexing: 
+    output_TFLAG!i_dim = input_TFLAG!i_dim
+  end do
+  copy_VarAtts(input_TFLAG, output_TFLAG)
+
+;   ; start debugging----------------------------------------------------
+;   print(this_fn+": pre-population: printVarSummary(output_TFLAG)==")
+;   printVarSummary(output_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+  ;;; iterate dims=(TSTEP, VAR, DATE-TIME)
+;   do i_tstep = 0, 0
+  do i_tstep = 0, output_TFLAG_tsteps_n-1
+;     do i_var = 0, 0 ; is all we're gonna do, but ...
+    do i_var = 0, output_TFLAG_vars_n-1
+      ; set DATE-TIME structure: see https://github.com/TomRoche/ioapi "For example"
+      if (i_tstep .eq. 24) then
+        output_TFLAG(i_tstep, i_var, :) = (/ output_TFLAG_day +1, 0 /)
+      else
+        output_TFLAG(i_tstep, i_var, :) = (/ output_TFLAG_day, i_tstep * 10000 /)
+      end if
+    end do ; i_var = 0, input_TFLAG_vars_n-1
+  end do ; i_tstep = 0, input_TFLAG_tsteps_n-1
+
+;   ; start debugging----------------------------------------------------
+;   print(this_fn+": post-population: print(output_TFLAG)==")
+;   print(output_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+;----------------------------------------------------------------------
+; write hourly output
+;----------------------------------------------------------------------
+;;; see procedure outline @ http://www.ncl.ucar.edu/Applications/method_2.shtml
+
+  ;;; 1. Get output file handle.
+  output_fh = addfile(output_fp, "c")
+
+  ;; 1.1. (optional) declare file definition mode
+  setfileoption(output_fh, "DefineMode", True)
+
+  ;;; 2. Define global/file attributes.
+
+  ;; 2.1. Copy global/file attributes from input: see Example 3 in http://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml
+  input_fa_names = getvaratts(input_fh)
+  if (.not. all(ismissing(input_fa_names))) then
+    do i=0, dimsizes(input_fa_names)-1
+      fa_name = input_fa_names(i)
+; should all be the same, possibly excepting :TSTEP?
+;       if (fa_name .eq. "NVARS") then
+;         output_fh@$fa_name$ = 1 ; |data variables|
+;       else if (fa_name .eq. "VAR-LIST") then
+;         ; IOAPI-style
+;         output_fh@$fa_name$ = pad_blanks_to_length(output_dv_name, 16)
+;       else
+        output_fh@$fa_name$ = input_fh@$fa_name$
+;       end if
+;       end if
+;       end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+    end do
+  end if
+
+  ;; 2.2. Overwrite our own global/file attributes.
+
+  ;;; 3. Define {dimensions, coordinate variables} and their dimensionality
+
+  ;;; don't copy all from input? we only want those needed for output_arr?
+
+  ;;; 3.1. Define output dimensions.
+  ;; names from input, sizes from output
+  output_fh_dims_names = getvardims(input_fh)
+  output_fh_dims_n = dimsizes(output_fh_dims_names)
+
+  output_fh_dims_sizes = new(output_fh_dims_n, integer, "No_FillValue")
+  do i_dim = 0 , output_fh_dims_n-1
+    i_dim_name = output_fh_dims_names(i_dim)
+    if      (i_dim_name .eq. "COL") then
+      output_fh_dims_sizes(i_dim) = output_dv_cols_n
+    else if (i_dim_name .eq. "DATE-TIME") then
+      output_fh_dims_sizes(i_dim) = output_TFLAG_datetimes_n
+    else if (i_dim_name .eq. "LAY") then
+      output_fh_dims_sizes(i_dim) = output_dv_layers_n
+    else if (i_dim_name .eq. "ROW") then
+      output_fh_dims_sizes(i_dim) = output_dv_rows_n
+    else if (i_dim_name .eq. "TSTEP") then
+      output_fh_dims_sizes(i_dim) = output_TFLAG_tsteps_n
+    else if (i_dim_name .eq. "VAR") then
+      output_fh_dims_sizes(i_dim) = output_TFLAG_vars_n
+    end if  
+    end if  
+    end if  
+    end if  
+    end if  
+    end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+  end do
+
+  ;; Should any dimensions be unlimited? Here, no.
+  output_fh_dims_unlim = new(output_fh_dims_n, logical, "No_FillValue")
+  do i=0 , output_fh_dims_n-1
+    output_fh_dims_unlim(i) = False
+  end do
+
+  filedimdef(output_fh,         \
+    (/ output_fh_dims_names /), \
+    (/ output_fh_dims_sizes /), \
+    (/ output_fh_dims_unlim /)  \
+  )
+
+; Skip this for now: there are no coordvars in the input, and we're more focused on IOAPI/VERDI compliance
+
+;   ;;; 3.2. Create output coordinate variables (if they exist in the input).
+;   ;;; procedure adapted from contributed.ncl::copy_VarCoords
+;   if (.not. all(ismissing(input_dv_dims_names))) then
+; ; Skip this for our case: we explicitly want |output dimensions| < |input dimensions|
+; ;     if (all(input_dv_dims_sizes(:output_dv_dims_n-1) .eq. output_dv_dims_sizes)) then
+;     do i=0 /), output_dv_dims_n-1
+; ; iterate over the output dimension names, which omit the first input dimension (name)
+;       if (.not. ismissing(output_dv_dims_names(i))) then
+;         output_arr!i = input_dv!(i+1)
+;         if (iscoord(input_dv, input_dv!(i+1)))
+;           output_arr&$output_arr!i$ = input_dv&$input_dv!i$
+;         end if
+;       end if
+;     end do
+; ;   else
+; ;     print(this_fn+": ERROR: dimension sizes do not match: input="+input_dv_dims_sizes+", output="+output_dv_dims_sizes)
+; ;   end if ; checking in and out dim sizes
+;   end if ; .not. all(ismissing(input_dv_dims_names))
+
+  ;;; 4. Setup output variable(s)
+
+  ;; 4.1. Define dimensionality of the output variable(s)
+  ; for VERDI/IOAPI-compatibility, do fake "datavar"=TFLAG
+  ; Try to make TFLAG the first "datavar"? Not sure it matters, but in the IOAPI files I've `ncdump`ed, it always is.
+  ; define order appears to determine `ncdump -h` order
+
+  filevardef(output_fh, output_TFLAG_name, output_TFLAG_type, output_TFLAG_dims_names)
+  ; now do the real data
+  filevardef(output_fh, output_dv_name, output_dv_type, output_dv_dims_names)
+
+  ;; 4.2. Define output variable(s) attributes
+
+  ; 4.2.1. Copy output variable(s) attributes
+  filevarattdef(output_fh, output_TFLAG_name, input_TFLAG)
+  filevarattdef(output_fh, output_dv_name, input_dv)
+
+  ;; 4.2.2. Overwrite our "real" output datavar(s) attributes
+;   output_fh->$output_dv_name$@units = output_dv_units
+;   output_fh->$output_dv_name$@long_name = output_dv_long_name
+;   output_fh->$output_dv_name$@var_desc = output_dv_var_desc
+  
+  ;; 4.3. (optional) Exit file definition mode
+  setfileoption(output_fh, "DefineMode", False)
+
+  ;;; 5. Write data to output variable(s)
+  output_fh->$output_TFLAG_name$ = (/ output_TFLAG /)
+  output_fh->$output_dv_name$ = (/ output_arr /)
+
+  ;;; 6. Close output filehandle to flush file: see http://www.ncl.ucar.edu/Support/talk_archives/2012/2196.html
+  delete(output_fh)
+
+; get/run VERDI command in driver rather than here, because driver moves the output before the process we launch here can get to it
+; ;----------------------------------------------------------------------
+; ; verify VERDI compatibility? "pure" output
+; ;----------------------------------------------------------------------
+; 
+;   ;;; TODO: integrate VERDI into bash_utilities.sh
+;   ;;; plot our datavar of interest with VERDI
+;   output_dv_name_VERDI = output_dv_name+"[1]"
+;   ;; VERDI really wants an FQP
+;   output_fp_VERDI = systemfunc("readlink -f "+output_fp)
+;   output_dv_VERDI_cmd = "verdi -f "+output_fp_VERDI+" -s "+output_dv_name_VERDI+" -gtype tile &"
+; 
+;   print(str_get_nl() + \
+;     this_fn+"; debugging IOAPI: about to run '"+output_dv_VERDI_cmd+"'" + str_get_nl())
+;   system(output_dv_VERDI_cmd)
+; ;   print(str_get_nl() + \
+; ;     this_fn+"; not debugging IOAPI: skipping '"+output_dv_VERDI_cmd+"'" + str_get_nl())
+
+end ; make_hourlies.ncl

File uber_driver.sh

            ./combine_EDGAR_and_EPIC_emissions.r \
            ./compute_EPIC_emissions.r \
            ./distill_EPIC_emittivities.ncl \
+           ./make_hourlies.ncl \
            ./read_BELD_CSV.r \
            ./regrid_vis.r \
            ./reunit_vis_EDGAR_global.r "