Commits

Tom Roche  committed f6720fe

all EPIC-related emissions VERDIfied

Upgraded combine_EDGAR_and_EPIC_emissions.r to make its output IOAPI-compliant sufficient to launch VERDI on it from the script. Thus now all of

* the "purified" EPIC emittivities https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/5yravg_20111219_pure.nc.gz

* the EPIC emissions (emittivities * BELD crop coverages) https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/5yravg_20111219_pure_emis.nc

* the combined (EPIC+EDGAR, type4) emissions https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/5yravg_20111219_pure_emis_both.nc

Also:

* for greater ease of debugging and testing, both {combine_EDGAR_and_EPIC_emissions.r, compute_EPIC_emissions.r} load environment variables (R::Sys.getenv) with 'unset=NA'

* added function=get.ncdf4.prec to compute_EPIC_emissions.r to workaround apparent bug with types/precisions in ncdf4 (TODO: move to regrid_utils)
* homogenized strings in AQMEII_driver.sh, causing need to tweak regrid_vis.r

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

TODO:
* handle {name, path to} VERDI in regrid_utils/bash_utilities.sh
* move R functions from {combine_EDGAR_and_EPIC_emissions.r, compute_EPIC_emissions.r} to regrid_utils/R_utilities.r
* 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
* AQMEII_driver.sh still has much uncalled cruft from ancestors=EDGAR_driver.sh, GFED_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 6d55d0c

Comments (0)

Files changed (4)

File AQMEII_driver.sh

 # IOAPI pads varattr=long_name to length=16 with trailing spaces
 AQMEIINA_TFLAG_LONG_NAME="$(printf '%-16s' ${AQMEIINA_TFLAG_NAME})"
 # IOAPI pads varattr=units to length=16 with trailing spaces
-AQMEIINA_DV_UNITS="$(printf '%-16s' '<YYYYDDD,HHMMSS>')"
+AQMEIINA_TFLAG_UNITS="$(printf '%-16s' '<YYYYDDD,HHMMSS>')"
 # IOAPI pads varattr=var_desc to length=80 with trailing spaces
-AQMEIINA_DV_VAR_DESC="$(printf '%-80s' 'Timestep-valid flags:  (1) YYYYDDD or (2) HHMMSS')"
+AQMEIINA_TFLAG_VAR_DESC="$(printf '%-80s' 'Timestep-valid flags:  (1) YYYYDDD or (2) HHMMSS')"
 
 AQMEIINA_DIM_LAYER_NAME='LAY'
 AQMEIINA_DIM_LAYER_LONG_NAME='index of layers above surface'
 export AQMEII_EPIC_N2O_EMITT_DATAVAR_VAR_DESC="$(printf '%-80s' 'N2O emittivities')"
 
 ## "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
-export AQMEII_EPIC_N2O_EMITT_TFLAG_NAME='TFLAG'
+export AQMEII_EPIC_N2O_EMITT_TFLAG_NAME="${AQMEIINA_TFLAG_NAME}"
 # with its own special dimensions
-export AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_DATETIME_NAME='DATE-TIME'
-export AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_VAR_NAME='VAR'
+export AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_DATETIME_NAME="${AQMEIINA_DIM_DATETIME_NAME}"
+export AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_VAR_NAME="${AQMEIINA_DIM_VAR_NAME}"
 
 # get these from driver
 # input.epic.TFLAG.units <- input.epic.TFLAG$units
 export AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_LONG_NAME="${AQMEIINA_DV_LONG_NAME}"
 export AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_UNITS="${AQMEIINA_DV_UNITS}"
 export AQMEII_EPIC_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_EPIC_N2O_EMISS_TFLAG_NAME="${AQMEIINA_TFLAG_NAME}"
 export AQMEII_EPIC_N2O_EMISS_TFLAG_ATTR_LONG_NAME="${AQMEIINA_TFLAG_LONG_NAME}"
 export AQMEII_EPIC_N2O_EMISS_TFLAG_ATTR_UNITS="${AQMEIINA_TFLAG_UNITS}"
 export AQMEII_EPIC_N2O_EMISS_TFLAG_ATTR_VAR_DESC="${AQMEIINA_TFLAG_VAR_DESC}"
 
 ## file/global attributes
-# don't need this ...
-# export AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS='1' # output semantics: "layer" == vertical, not crop
-# but *do* need this ...
+export AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS='1' # output semantics: "layer" == vertical, not crop
 export AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS_TYPE='integer' # 'short' -> `:NLAYS = 1s ;`
 export AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST="${AQMEIINA_DV_LONG_NAME}"
 export AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST_TYPE='character'
 export AQMEII_BOTH_N2O_EMISS_DATAVAR_ATTR_LONG_NAME="${AQMEIINA_DV_LONG_NAME}"
 export AQMEII_BOTH_N2O_EMISS_DATAVAR_ATTR_UNITS="${AQMEIINA_DV_UNITS}"
 export AQMEII_BOTH_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_BOTH_N2O_EMISS_TFLAG_NAME="${AQMEIINA_TFLAG_NAME}"
+export AQMEII_BOTH_N2O_EMISS_TFLAG_ATTR_LONG_NAME="${AQMEIINA_TFLAG_LONG_NAME}"
+export AQMEII_BOTH_N2O_EMISS_TFLAG_ATTR_UNITS="${AQMEIINA_TFLAG_UNITS}"
+export AQMEII_BOTH_N2O_EMISS_TFLAG_ATTR_VAR_DESC="${AQMEIINA_TFLAG_VAR_DESC}"
 
 ## file/global attributes
 # export AQMEII_BOTH_N2O_EMISS_ATTR_NLAYS='1' # output semantics: "layer" == vertical, not crop
 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_TYPE='text' # or 'character': ncdf4 types
+export AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC_TYPE='character'
 
 ## visualization
 AQMEII_BOTH_N2O_EMISS_PDF_FN="${AQMEII_BOTH_N2O_EMISS_ROOT}.pdf"
   'BELD_CSV_to_RDS' \
   'distill_EPIC' \
   'compute_EPIC' \
+  'combine_emis' \
   'teardown' \
 ; do
   echo -e "\n$ ${THIS_FN}::main loop::${CMD}\n"
     * [raster][]
     * [rasterVis][]
 
+Visualizations of datasets are provided
+
+* statistically (summary statistics of datasets are generated and presented in console)
+* by plots (saved as PDF, displayed by launching the user-configured PDF viewer)
+* by [VERDI][] (which also serves as a check on the [IOAPI][]-compatibility of the generated netCDF)
+
 [bash @ wikipedia]: http://en.wikipedia.org/wiki/Bash_%28Unix_shell%29
 [R @ wikipedia]: http://en.wikipedia.org/wiki/R_%28programming_language%29
 [NCL @ wikipedia]: http://en.wikipedia.org/wiki/NCAR_Command_Language
 [EDGAR 4.2 overview]: http://edgar.jrc.ec.europa.eu/overview.php?v=42
 [AQMEII-NA]: https://github.com/TomRoche/cornbeltN2O/wiki/AQMEII-North-American-domain
 [BELD overview]: TODO!
-[EPIC overview]: TODO!
+[EPIC overview]: https://github.com/TomRoche/cornbeltN2O/wiki/EPIC
 [epic_site_crops_0529_USA_2Ellen.csv.gz]: https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/epic_site_crops_0529_USA_2Ellen.csv.gz
 [R internals::Serialization Formats]: http://cran.r-project.org/doc/manuals/R-ints.html#Serialization-Formats
 [visualization of global EDGAR]: https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/v42_N2O_2008_IPCC_4C_4D.0.1x0.1_reunit.pdf
 [ncdf4]: http://cran.r-project.org/web/packages/ncdf4/
 [raster]: http://cran.r-project.org/web/packages/raster/
 [rasterVis]: http://cran.r-project.org/web/packages/rasterVis/
+[VERDI]: http://www.verdi-tool.org/
+[IOAPI]: http://www.baronams.com/products/ioapi/
 
 To run this code,
 

File combine_EDGAR_and_EPIC_emissions.r

 input_edgar_fp <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_FP', unset=NA) # path to output file
 
 # file/global attributes
-input_edgar_attr_nlays = as.numeric(Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_NLAYS'))
-input_edgar_attr_nlays_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_NLAYS_TYPE')
-input_edgar_attr_varlist = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_VARLIST')
-input_edgar_attr_varlist_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_VARLIST_TYPE')
-input_edgar_attr_filedesc = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_FILEDESC')
-input_edgar_attr_filedesc_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_FILEDESC_TYPE')
+input_edgar_attr_nlays = as.numeric(Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_NLAYS', unset=NA))
+input_edgar_attr_nlays_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_NLAYS_TYPE', unset=NA)
+input_edgar_attr_varlist = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_VARLIST', unset=NA)
+input_edgar_attr_varlist_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_VARLIST_TYPE', unset=NA)
+input_edgar_attr_filedesc = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_FILEDESC', unset=NA)
+input_edgar_attr_filedesc_type = Sys.getenv('AQMEII_EDGAR_N2O_REGRID_ATTR_FILEDESC_TYPE', unset=NA)
 
 # dimensions
 input_edgar_dim_layers_name <-  Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_LAYERS_NAME', unset=NA)
 input_edgar_dv_long_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DATAVAR_ATTR_LONG_NAME', unset=NA)
 input_edgar_dv_var_desc <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DATAVAR_ATTR_VAR_DESC', unset=NA)
 
-## EPIC N2O emissions over AQMEII-NA (actually CONUS)
+## EPIC N2O emissions over AQMEII-NA (actually CONUS): reunit-ed, monotonic, VERDI-compliant
 
 input_epic_fp <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_FP', unset=NA) # path to output file
 
 # file/global attributes
-input_epic_attr_nlays = as.numeric(Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS'))
-input_epic_attr_nlays_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS_TYPE')
-input_epic_attr_varlist = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST')
-input_epic_attr_varlist_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST_TYPE')
-input_epic_attr_filedesc = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC')
-input_epic_attr_filedesc_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC_TYPE')
+input_epic_attr_nlays = as.numeric(Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS', unset=NA))
+input_epic_attr_nlays_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS_TYPE', unset=NA)
+input_epic_attr_varlist = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST', unset=NA)
+input_epic_attr_varlist_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST_TYPE', unset=NA)
+input_epic_attr_filedesc = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC', unset=NA)
+input_epic_attr_filedesc_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC_TYPE', unset=NA)
 
 # dimensions
-input_epic_dim_layers_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_NAME', unset=NA)
-input_epic_dim_layers_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_ATTR_UNITS', unset=NA)
-input_epic_dim_layers_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_ATTR_LONG_NAME', unset=NA)
+input_epic_dim_datetimes_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_NAME', unset=NA)
+input_epic_dim_layers_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_NAME', unset=NA)
+input_epic_dim_layers_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_UNITS', unset=NA)
+input_epic_dim_layers_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_LONG_NAME', unset=NA)
+input_epic_dim_tsteps_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_NAME', unset=NA)
+input_epic_dim_vars_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_VAR_NAME', unset=NA)
 input_epic_dim_x_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_NAME', unset=NA)
 input_epic_dim_x_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_UNITS', unset=NA)
 input_epic_dim_x_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_LONG_NAME', unset=NA)
 input_epic_dv_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_UNITS', unset=NA)
 input_epic_dv_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_LONG_NAME', unset=NA)
 input_epic_dv_var_desc <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_VAR_DESC', unset=NA)
+# "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
+input_epic_TFLAG_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_TFLAG_NAME', unset=NA)
+# input_epic_TFLAG_var_desc <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_TFLAG_ATTR_VAR_DESC', unset=NA) # get/copy programmatically
 
 ## outputs
 
-## combined N2O emissions
+## combined N2O emissions, VERDI-compliant
 
 output_both_fp <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_FP', unset=NA) # path to output file
+
 # file/global attributes
-output_both_attr_nlays = as.numeric(Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_NLAYS'))
-output_both_attr_nlays_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_NLAYS_TYPE')
-output_both_attr_varlist = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST')
-output_both_attr_varlist_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST_TYPE')
-output_both_attr_filedesc = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC')
-output_both_attr_filedesc_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC_TYPE')
+output_both_attr_nlays = as.numeric(Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_NLAYS', unset=NA))
+output_both_attr_nlays_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_NLAYS_TYPE', unset=NA)
+output_both_attr_varlist = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST', unset=NA)
+output_both_attr_varlist_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_VARLIST_TYPE', unset=NA)
+output_both_attr_filedesc = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC', unset=NA)
+output_both_attr_filedesc_type = Sys.getenv('AQMEII_BOTH_N2O_EMISS_ATTR_FILEDESC_TYPE', unset=NA)
+
 # dimensions
+output_both_dim_datetimes_name <-  Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_DATETIME_NAME', unset=NA)
+output_both_dim_tsteps_name <-  Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_TSTEP_NAME', unset=NA)
+output_both_dim_vars_name <-  Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_VAR_NAME', unset=NA)
 output_both_dim_layers_name <-  Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_LAYERS_NAME', unset=NA)
 output_both_dim_layers_units <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_LAYERS_ATTR_UNITS', unset=NA)
 output_both_dim_layers_long_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_LAYERS_ATTR_LONG_NAME', unset=NA)
 output_both_dim_y_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_Y_NAME', unset=NA)
 output_both_dim_y_units <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_Y_ATTR_UNITS', unset=NA)
 output_both_dim_y_long_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DIM_Y_ATTR_LONG_NAME', unset=NA)
+
 # datavars
 output_both_dv_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DATAVAR_NAME', unset=NA)
 output_both_dv_units <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DATAVAR_ATTR_UNITS', unset=NA)
 output_both_dv_long_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DATAVAR_ATTR_LONG_NAME', unset=NA)
 output_both_dv_var_desc <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_DATAVAR_ATTR_VAR_DESC', unset=NA)
+# "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
+output_both_TFLAG_name <- Sys.getenv('AQMEII_BOTH_N2O_EMISS_TFLAG_NAME', unset=NA)
+# output_both_TFLAG_var_desc <-  Sys.getenv('AQMEII_BOTH_N2O_EMISS_TFLAG_ATTR_VAR_DESC', unset=NA) # get/copy programmatically
 
 # visualization
 
 # internal functions (defined)
 # ----------------------------------------------------------------------
 
+### TODO: move to regrid_utils: get.ncdf4.prec also used by compute_EPIC_emissions.r
+
+# workaround apparent bug in ncdf4
+get.ncdf4.prec <- function(
+  ncvar4 # object of type=ncvar4
+) {
+  # TODO: test argument
+  ncdf4.prec <- ncvar4$prec
+  if (ncdf4.prec == "int") {
+    return("integer")
+  } else {
+    return(ncdf4.prec)
+  }
+}
+
 ### Test if argument is zero within tolerance
 is.zero.tol <- function (val, tol = getOption("tolerance")) {
   if (is.null(tol)) {
 
 library(ncdf4)
 
+### EDGAR N2O emissions
+
+input.edgar.fh <- ncdf4::nc_open(input_edgar_fp, write=FALSE, readunlim=FALSE)
+input.edgar.arr <- ncdf4::ncvar_get(input.edgar.fh, varid=input_edgar_dv_name)
+input.edgar.dv <- input.edgar.fh$var[[input_edgar_dv_name]] # note double brackets!
+input.edgar.dv.type <- input.edgar.dv$prec
+input.edgar.dv.mv <- input.edgar.dv$missval
+input.edgar.dim.x <- input.edgar.fh$dim[[input_edgar_dim_x_name]]
+input.edgar.dim.y <- input.edgar.fh$dim[[input_edgar_dim_y_name]]
+
 ### EPIC N2O emissions
 
 input.epic.fh <- ncdf4::nc_open(input_epic_fp, write=FALSE, readunlim=FALSE)
 # input.epic.dv <- ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name)
 # `ncdf4::ncvar_get` gets data, not metadata
-input.epic.arr <- ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name)
+input.epic.dv.arr <- ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name)
 input.epic.dv <- input.epic.fh$var[[input_epic_dv_name]] # note double brackets!
 
 ## get input datavar type, missing value, long_name, units useful for output
 # input.epic.dv.long_name <- input.epic.dv$long_name
 # input.epic.dv.var_desc <- input.epic.dv$var_desc
 
-## get input dimensions useful for output (i.e., of type=ncdf4, not type=list)
+## "fake" datavar=TFLAG, for IOAPI and (IIUC) VERDI
+input.epic.TFLAG.arr <- ncdf4::ncvar_get(input.epic.fh, varid=input_epic_TFLAG_name)
+input.epic.TFLAG.dims <- dim(input.epic.TFLAG.arr)
+input.epic.TFLAG.dims.n <- length(input.epic.TFLAG.dims)
+input.epic.TFLAG <- input.epic.fh$var[[input_epic_TFLAG_name]]
+input.epic.TFLAG.type <- input.epic.TFLAG$prec
+input.epic.TFLAG.mv <- input.epic.TFLAG$missval
+
+### get input dimensions useful for output (i.e., of type=ncdf4, not type=list)
 # # print(input.epic.fh$dim)
 # output.epic.dim.y <- input.epic.fh$dim[2]
 # output.epic.dim.x <- input.epic.fh$dim[3]
-input.epic.dim.x <- input.epic.fh$dim[[input_epic_dim_x_name]]
 input.epic.dim.y <- input.epic.fh$dim[[input_epic_dim_y_name]]
-
-### EDGAR N2O emissions
-
-input.edgar.fh <- ncdf4::nc_open(input_edgar_fp, write=FALSE, readunlim=FALSE)
-input.edgar.arr <- ncdf4::ncvar_get(input.edgar.fh, varid=input_edgar_dv_name)
-input.edgar.dv <- input.edgar.fh$var[[input_edgar_dv_name]] # note double brackets!
-input.edgar.dv.type <- input.edgar.dv$prec
-input.edgar.dv.mv <- input.edgar.dv$missval
-input.edgar.dim.x <- input.edgar.fh$dim[[input_edgar_dim_x_name]]
-input.edgar.dim.y <- input.edgar.fh$dim[[input_edgar_dim_y_name]]
-
-### TODO: gotta set layers, otherwise output file gets written with
-### > :NLAYS = -2147483648 ;
-
-### TODO: gotta set timesteps (here and in ancestors :-(), otherwise output file gets written with
-### > :TSTEP = 240000 ;
+input.epic.dim.x <- input.epic.fh$dim[[input_epic_dim_x_name]]
+input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
+input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
+input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
+input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
 
 # ----------------------------------------------------------------------
 # get input and output dimensionalities
 ### missing_value may differ, however!
 ### and prefer metadata from EPIC (it's more IOAPI-ish)
 
-### Both EPIC input and output have same spatial grid, input also has layers=crops as leftmost dim (reversed by R)
-# input.epic.dv.dims <- dim(input.epic.dv)
-# > NULL
-input.edgar.dv.dims <- dim(input.edgar.arr)
-input.epic.dv.dims <- dim(input.epic.arr)
-base::stopifnot(input.edgar.dv.dims == input.epic.dv.dims)
-output.both.dv.dims <- input.edgar.dv.dims
+### EDGAR input
 
+input.edgar.dv.dims <- dim(input.edgar.arr)
 input.edgar.dv.dims.n <- length(input.edgar.dv.dims)
-input.epic.dv.dims.n <- length(input.epic.dv.dims)
-base::stopifnot(input.edgar.dv.dims.n == input.epic.dv.dims.n)
-output.both.dv.dims.n <- input.edgar.dv.dims.n
-
 input.edgar.dim.x.n <- input.edgar.dim.x$len
-input.epic.dim.x.n <- input.epic.dim.x$len
-base::stopifnot(input.edgar.dim.x.n == input.epic.dim.x.n)
-output.both.dim.x.n <- input.edgar.dim.x.n
-
 input.edgar.dim.y.n <- input.edgar.dim.y$len
+input.edgar.dv.cells.n <- input.edgar.dim.y.n * input.edgar.dim.x.n # spatial only
+
+### Both EPIC input and output have same spatial grid, input also has layers=1 and tsteps (probably=1)
+
+### EPIC input
+
+## EPIC input dims
+
+input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
+input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
+input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
+input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
+input.epic.dim.x <- input.epic.fh$dim[[input_epic_dim_x_name]]
+input.epic.dim.y <- input.epic.fh$dim[[input_epic_dim_y_name]]
+
+input.epic.dim.datetimes.n <- input.epic.dim.datetimes$len
+input.epic.dim.layers.n <- input.epic.dim.layers$len
+input.epic.dim.tsteps.n <- input.epic.dim.tsteps$len
+input.epic.dim.vars.n <- input.epic.dim.vars$len
+input.epic.dim.x.n <- input.epic.dim.x$len
 input.epic.dim.y.n <- input.epic.dim.y$len
-base::stopifnot(input.edgar.dim.y.n == input.epic.dim.y.n)
-output.both.dim.y.n <- input.edgar.dim.y.n
 
-input.edgar.dv.cells.n <- input.edgar.dim.y.n * input.edgar.dim.x.n # spatial only
-input.epic.dv.cells.n <- input.edgar.dv.cells.n
-output.both.dv.cells.n <- input.edgar.dv.cells.n
+## EPIC input "real" datavar
+
+input.epic.dv.dims <- dim(input.epic.dv.arr)
+input.epic.dv.dims.n <- length(input.epic.dv.dims)
+# spatial cells only: ignore layers, tsteps
+input.epic.dv.cells.n <- input.epic.dim.y.n * input.epic.dim.x.n
+
+## EPIC input TFLAG
+
+input.epic.TFLAG.dims <- dim(input.epic.TFLAG.arr)
+input.epic.TFLAG.dims.n <- length(input.epic.TFLAG.dims)
+
+### combined output
+
+output.both.dim.datetimes.n <- input.epic.dim.datetimes.n
+output.both.dim.layers.n <- input.epic.dim.layers.n
+output.both.dim.tsteps.n <- input.epic.dim.tsteps.n
+output.both.dim.vars.n <- input.epic.dim.vars.n
+output.both.dim.x.n <- input.epic.dim.x.n
+output.both.dim.y.n <- input.epic.dim.y.n
+
+output.both.dv.cells.n <- input.epic.dv.cells.n
+output.both.dv.dims <-
+  c(output.both.dim.x.n, output.both.dim.y.n, output.both.dim.layers.n, output.both.dim.tsteps.n)
+output.both.dv.dims.n <- length(output.both.dv.dims)
 
 # ----------------------------------------------------------------------
 # read inputs
 
 ### EPIC input
 
-input.epic.read.start <- rep(1, input.epic.dv.dims.n)   # start=(1,1,1,...)
-input.epic.read.count <- input.epic.dv.dims             # count=(n,n,n,...)
-input.epic.arr <-
-  ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name, start=input.epic.read.start, count=input.epic.read.count)
+## TODO: "Pierce-style read" (see help(ncvar_get)#Examples): protects against data overload (somewhat)
+## but more complicated to implement. Copy code from (github)ioapi-hack-R/computeCropSum.r
+## Normally one reads one timestep at a time, but EPIC input has no timesteps--
+## other than the implicit/degenerate model year ... which, for VERDI-compliance, it appears we must use :-(
+## So gotta deal with implicit/degenerate IOAPI dimensions ...
+
+# so gotta handle both LAY=1 and TSTEP=1: remembering timelike dim is always the LAST dimension!
+if      (input.epic.dv.dims.n < 2) {
+  # TODO: throw
+  stop(sprintf('%s: ERROR: dimensionality=%i < 2 of EPIC input datavar=%s, dimensions are',
+    this_fn, input.epic.dv.dims.n, input_epic_dv_name))
+  print(input.epic.dv.dims)
+
+} else if (input.epic.dv.dims.n == 2) {
+  if ((input.epic.dim.tsteps.n != 1) && (input.epic.dim.lay.n != 1)) {
+    stop(sprintf('%s: ERROR: dimensionality=1 of EPIC input datavar=%s, but dim(TSTEP)=%i and dim(LAY)=%i',
+      this_fn, input_epic_TFLAG_name, input.epic.dim.tsteps.n, input.epic.dim.lays.n))
+  }
+  # append TSTEP, LAY -> (TSTEP, LAY, ROW, COL) ... but reversed :-(
+  input.epic.dv.dims <- c(input.epic.dv.dims, 1, 1)
+  input.epic.dv.dims.n <- 4
+  input.epic.dv.read.start <- rep(1, input.epic.dv.dims.n)   # start=(1,1,1,...)
+  input.epic.dv.read.count <- input.epic.dv.dims             # count=(n,n,n,...)
+
+# start debugging-------------------------------------------
+  cat(sprintf('\n%s: dimensionality=%i (was 2) of EPIC input datavar=%s, dimensions are',
+    this_fn, input.epic.dv.dims.n, input_epic_dv_name))
+  print(input.epic.dv.dims)
+#   end debugging-------------------------------------------
+
+} else if (input.epic.dv.dims.n == 3) {
+  if        (input.epic.dim.tsteps.n == 1) {
+    # append TSTEPs
+    input.epic.dv.dims <- c(input.epic.dv.dims, 1)
+  } else if  (input.epic.dim.lays.n == 1) {
+    # insert LAYs
+    input.epic.dv.dims <-
+      c(input.epic.dim.x.n, input.epic.dim.y.n, 1, input.epic.dim.tsteps.n)
+  } else {
+    stop(sprintf('%s: ERROR: dimensionality=3 of EPIC input datavar=%s, but dim(TSTEP)==%i and dim(LAY)==%i',
+      this_fn, input_epic_dv_name, input.epic.dim.tsteps.n, input.epic.dim.lays.n))
+  } # end testing LAYs, TSTEPs
+  input.epic.dv.dims.n <- 4
+  input.epic.dv.read.start <- rep(1, input.epic.dv.dims.n)
+  input.epic.dv.read.count <- input.epic.dv.dims
+
+# start debugging-------------------------------------------
+  cat(sprintf('\n%s: dimensionality=%i (was 3) of EPIC input datavar=%s, dimensions are',
+    this_fn, input.epic.dv.dims.n, input_epic_dv_name))
+  print(input.epic.dv.dims)
+#   end debugging-------------------------------------------
+
+} else if (input.epic.dv.dims.n == 4) {
+# start debugging-------------------------------------------
+  cat(sprintf('\n%s: dimensionality=4 of EPIC input datavar=%s, dimensions are',
+    this_fn, input_epic_dv_name))
+  print(input.epic.dv.dims)
+#   end debugging-------------------------------------------
+  input.epic.dv.read.start <- rep(1, input.epic.dv.dims.n)
+  input.epic.dv.read.count <- input.epic.dv.dims
+
+} else {
+  # TODO: throw
+  stop(sprintf('%s: ERROR: dimensionality=%i > 4 of EPIC input datavar=%s, dimensions are',
+    this_fn, input.epic.dv.dims.n, input_epic_dv_name))
+  print(input.epic.dv.dims)
+} # end testing input.epic.dv.dims.n
+
+### Get input data
+input.epic.dv.arr <-
+  ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name, start=input.epic.dv.read.start, count=input.epic.dv.read.count)
+## Once again, caught between IOAPI (which demands tstep=1) and netCDF (which strips tstep=1) :-(
+dim(input.epic.dv.arr) <- input.epic.dv.dims
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(input.epic.arr)==', this_fn))
-print(dim(input.epic.arr))
+base::cat(base::sprintf('\n%s: dim(input.epic.dv.arr)==', this_fn))
+print(dim(input.epic.dv.arr))
 stats.to.stdout(
-  data=input.epic.arr,
+  data=input.epic.dv.arr,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: EPIC input emissions=', this_fn)
 )
 
 # combine_EDGAR_and_EPIC_emissions.r: EPIC input emissions=
-#       cells=137241
-#       obs=40764
-#       min=0
-#       q1=1.02e-05
-#       med=0.000171
-#       mean=0.00138
-#       q3=0.00111
-#       max=0.283
-#       sum=56.2
+#         cells=137241
+#         obs=40763
+#         min=0
+#         q1=1.02e-05
+#         med=0.000171
+#         mean=0.00137
+#         q3=0.00111
+#         max=0.0757
+#         sum=56         [after removing outlier]
 #   end debugging-----------------------------------------------------
 
-### how does EPIC compare to EDGAR?
+### input TFLAG
 
-epic.gt.edgar.arr <- input.epic.arr - input.edgar.arr
+## ... gotta deal with implicit/degenerate TSTEP==1 again!
+## Remembering timelike dim is always the LAST dimension!
+if      (input.epic.TFLAG.dims.n < 1) {
+  # TODO: throw
+  stop(sprintf('%s: ERROR: dimensionality=%i < 1 of EPIC input TFLAG=%s',
+    this_fn, input.epic.TFLAG.dims.n, input_epic_TFLAG_name))
 
-# start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(epic.gt.edgar.arr)==', this_fn))
-print(dim(epic.gt.edgar.arr))
-stats.to.stdout(
-  data=epic.gt.edgar.arr,
-  sig.digs=sigdigs,
-  title=base::sprintf('\n%s: EPIC - EDGAR=', this_fn)
-)
+} else if (input.epic.TFLAG.dims.n == 1) {
+  if ((input.epic.dim.tsteps.n != 1) && (input.epic.dim.var.n != 1)) {
+    stop(sprintf('%s: ERROR: dimensionality=1 of EPIC input TFLAG=%s, but dim(TSTEP)=%i and dim(VAR)=%i',
+      this_fn, input_epic_TFLAG_name, input.epic.dim.tsteps.n, input.epic.dim.vars.n))
+  }
+  # append TSTEP, VAR -> (TSTEP, VAR, DATE-TIME) ... but reversed :-(
+  input.epic.TFLAG.dims <- c(input.epic.TFLAG.dims, 1, 1)
+  input.epic.TFLAG.dims.n <- 3
+# start debugging-------------------------------------------
+  cat(sprintf('\n%s: dimensionality=%i (was 1) of EPIC input TFLAG=%s, dimensions are',
+    this_fn, input.epic.TFLAG.dims.n, input_epic_TFLAG_name))
+  print(input.epic.TFLAG.dims)
+#   end debugging-------------------------------------------
+
+} else if (input.epic.TFLAG.dims.n == 2) {
+  # which |dim|==1 ???
+  if        (input.epic.dim.tsteps.n == 1) {
+    input.epic.TFLAG.dims <-
+      c(input.epic.dim.datetime.n, input.epic.dim.vars.n, 1)
+  } else if (input.epic.dim.var.n == 1) {
+    input.epic.TFLAG.dims <-
+      c(input.epic.dim.datetime.n, 1, input.epic.dim.tsteps.n)
+  }
+  input.epic.TFLAG.dims.n <- 3
+# start debugging-------------------------------------------
+  cat(sprintf('\n%s: dimensionality=%i (was 2) of EPIC input TFLAG=%s, dimensions are',
+    this_fn, input.epic.TFLAG.dims.n, input_epic_TFLAG_name))
+  print(input.epic.TFLAG.dims)
+#   end debugging-------------------------------------------
+
+} else if (input.epic.TFLAG.dims.n == 3) {
+
+  # life is good ...
+  cat(sprintf('\n%s: dimensionality=3 of EPIC input TFLAG=%s, dimensions are',
+    this_fn, input_epic_TFLAG_name))
+  print(input.epic.TFLAG.dims)
+
+} else {
+  # TODO: throw
+  stop(sprintf('%s: ERROR: dimensionality=%i > 3 of EPIC input TFLAG=%s',
+    this_fn, input.epic.TFLAG.dims.n, input_epic_TFLAG_name))
+} # end testing input.epic.TFLAG.dims.n
+input.epic.TFLAG.read.start <- rep(1, input.epic.TFLAG.dims.n)
+input.epic.TFLAG.read.count <- input.epic.TFLAG.dims
+
+# # ----------------------------------------------------------------------
+# # compare EPIC to EDGAR
+# # ----------------------------------------------------------------------
+
+# ### TODO: make dimensions of EPIC input match EDGAR input
+# epic.gt.edgar.arr <- input.epic.dv.arr - input.edgar.arr
 
-# combine_EDGAR_and_EPIC_emissions.r: EPIC - EDGAR=
-#         cells=137241
-#         obs=40764
-#         min=-0.0856
-#         q1=-0.00539
-#         med=-0.00229
-#         mean=-0.00249
-#         q3=-6.25e-06
-#         max=0.28
-#         sum=-102
-#   end debugging-----------------------------------------------------
+# # start debugging-----------------------------------------------------
+# base::cat(base::sprintf('\n%s: dim(epic.gt.edgar.arr)==', this_fn))
+# print(dim(epic.gt.edgar.arr))
+# stats.to.stdout(
+#   data=epic.gt.edgar.arr,
+#   sig.digs=sigdigs,
+#   title=base::sprintf('\n%s: EPIC - EDGAR=', this_fn)
+# )
+
+# # combine_EDGAR_and_EPIC_emissions.r: EPIC - EDGAR=
+# #         cells=137241
+# #         obs=40764
+# #         min=-0.0856
+# #         q1=-0.00539
+# #         med=-0.00229
+# #         mean=-0.00249
+# #         q3=-6.25e-06
+# #         max=0.28
+# #         sum=-102
+# #   end debugging-----------------------------------------------------
+
+# ### TODO: *plot* differences between EPIC input match EDGAR input
 
 # ----------------------------------------------------------------------
 # setup output data container(s)
 # ----------------------------------------------------------------------
 
-output.both.arr <- array(NA, output.both.dv.dims)
-output.both.write.start <- rep(1, output.both.dv.dims.n)   # start=(1,1,1,...)
-output.both.write.count <- output.both.dv.dims             # count=(n,n,n,...)
+### TFLAG
+
+output.both.TFLAG.dims <- input.epic.TFLAG.dims # (TSTEP, VAR, DATE-TIME), except ...
+output.both.TFLAG.dims.n <- length(output.both.TFLAG.dims)
+output.both.TFLAG.arr <- input.epic.TFLAG.arr
+output.both.TFLAG.write.start <- input.epic.TFLAG.read.start
+output.both.TFLAG.write.count <- input.epic.TFLAG.read.count
+
+# # start debugging-------------------------------------------
+# cat(sprintf('%s: output.both.TFLAG.dims==', this_fn))
+# print(output.both.TFLAG.dims)
+# cat(sprintf('%s: output.both.TFLAG.write.start==', this_fn))
+# print(output.both.TFLAG.write.start)
+# cat(sprintf('%s: output.both.TFLAG.write.count==', this_fn))
+# print(output.both.TFLAG.write.count)
+# #   end debugging-------------------------------------------
+
+### the real datavar
+
+output.both.dv.arr <- array(NA, output.both.dv.dims)
+output.both.dv.write.start <- rep(1, output.both.dv.dims.n)   # start=(1,1,1,...)
+output.both.dv.write.count <- output.both.dv.dims             # count=(n,n,n,...)
 
 # # start debugging-----------------------------------------------------
 # base::cat(base::sprintf('\n%s: output.both.write.start==', this_fn))
 base::cat(base::sprintf('\n%s: about to combine EDGAR and EPIC crop emissions\n', this_fn))
 #   end debugging-----------------------------------------------------
 
-### type1: use values from EDGAR only where EPIC value is missing/NA
+## regarding copy criteria, see https://stat.ethz.ch/pipermail/r-help/2013-May/353195.html
+## much more elegant than (all-the-way) nested `for` loops. thanks, Bill Dunlap!
+
+### ASSERT: EDGAR is a simple spatial grid (ROW, COL), while EPIC is a grid of the same spatiality possibly embedded in LAY and TSTEP.
 
-## newstyle! see https://stat.ethz.ch/pipermail/r-help/2013-May/353195.html
-## much more elegant than nested `for` loops. thanks, Bill Dunlap!
+### type1: use values from EDGAR only where EPIC value is missing/NA
 
-output.both.arr.type1 <- input.epic.arr
-# (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
-copy.criterion.type1 <- is.na(input.epic.arr) & !is.na(input.edgar.arr)
-output.both.arr.type1[copy.criterion.type1] <- input.edgar.arr[copy.criterion.type1]
+output.both.dv.arr.type1 <- input.epic.dv.arr
+for (i.tstep in 1:output.both.dim.tsteps.n) {
+  for (i.lay in 1:output.both.dim.layers.n) {
+    input.epic.grid <- input.epic.dv.arr[,,i.lay,i.tstep]
+    output.grid <- input.epic.grid
+    # (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
+    copy.criterion.type1 <-
+      is.na(input.epic.grid) & !is.na(input.edgar.arr)
+    output.grid[copy.criterion.type1] <- input.edgar.arr[copy.criterion.type1]
+    output.both.dv.arr.type1[,,i.lay,i.tstep] <- output.grid
+  }
+}
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(output.both.arr.type1)==', this_fn))
-print(dim(output.both.arr.type1))
+base::cat(base::sprintf('\n%s: dim(output.both.dv.arr.type1)==', this_fn))
+print(dim(output.both.dv.arr.type1))
 stats.to.stdout(
-  data=output.both.arr.type1,
+  data=output.both.dv.arr.type1,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: combined (type1) output emissions=', this_fn)
 )
 #         med=3.34e-05
 #         mean=0.00115
 #         q3=0.00109
-#         max=0.283
-#         sum=158
+#         max=0.0757
+#         sum=157
 
 # type1 < EDGAR: from above:
 # combine_EDGAR_and_EPIC_emissions.r: EDGAR input emissions=
 
 ### type2: use values from EDGAR only where EPIC value is zero or missing/NA (and EDGAR is not)
 
-output.both.arr.type2 <- input.epic.arr
-# (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
-copy.criterion.type2 <- is.arr.na.or.zero(input.epic.arr) & !is.arr.na.or.zero(input.edgar.arr)
-output.both.arr.type2[copy.criterion.type2] <- input.edgar.arr[copy.criterion.type2]
+output.both.dv.arr.type2 <- input.epic.dv.arr
+for (i.tstep in 1:output.both.dim.tsteps.n) {
+  for (i.lay in 1:output.both.dim.layers.n) {
+    input.epic.grid <- input.epic.dv.arr[,,i.lay,i.tstep]
+    output.grid <- input.epic.grid
+    # (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
+    copy.criterion.type2 <-
+      is.arr.na.or.zero(input.epic.grid) & !is.arr.na.or.zero(input.edgar.arr)
+    output.grid[copy.criterion.type2] <- input.edgar.arr[copy.criterion.type2]
+    output.both.dv.arr.type2[,,i.lay,i.tstep] <- output.grid
+  }
+}
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(output.both.arr.type2)==', this_fn))
-print(dim(output.both.arr.type2))
+base::cat(base::sprintf('\n%s: dim(output.both.dv.arr.type2)==', this_fn))
+print(dim(output.both.dv.arr.type2))
 stats.to.stdout(
-  data=output.both.arr.type2,
+  data=output.both.dv.arr.type2,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: combined (type2) output emissions=', this_fn)
 )
 #         min=0
 #         q1=3.14e-05
 #         med=0.00043
-#         mean=0.00171
+#         mean=0.0017
 #         q3=0.00237
-#         max=0.283
-#         sum=158
+#         max=0.0757
+#         sum=157
 
 # type2 < EDGAR: from above:
 # combine_EDGAR_and_EPIC_emissions.r: EDGAR input emissions=
 
 ### type3: use values from EPIC only where EDGAR value is zero or NA
 
-output.both.arr.type3 <- input.edgar.arr
-# (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
-copy.criterion.type3 <- is.arr.na.or.zero(input.edgar.arr) & !is.arr.na.or.zero(input.epic.arr)
-output.both.arr.type3[copy.criterion.type3] <- input.epic.arr[copy.criterion.type3]
+output.both.dv.arr.type3 <- input.epic.dv.arr # can't use EDGAR: dims don't match
+for (i.tstep in 1:output.both.dim.tsteps.n) {
+  for (i.lay in 1:output.both.dim.layers.n) {
+    input.epic.grid <- input.epic.dv.arr[,,i.lay,i.tstep]
+    output.grid <- input.edgar.arr
+    # (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
+    copy.criterion.type3 <-
+      is.arr.na.or.zero(input.edgar.arr) & !is.arr.na.or.zero(input.epic.grid)
+    output.grid[copy.criterion.type3] <- input.epic.grid[copy.criterion.type3]
+    output.both.dv.arr.type3[,,i.lay,i.tstep] <- output.grid
+  }
+}
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(output.both.arr.type3)==', this_fn))
-print(dim(output.both.arr.type3))
+base::cat(base::sprintf('\n%s: dim(output.both.dv.arr.type3)==', this_fn))
+print(dim(output.both.dv.arr.type3))
 stats.to.stdout(
-  data=output.both.arr.type3,
+  data=output.both.dv.arr.type3,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: combined (type3) output emissions=', this_fn)
 )
 # type3 vs type1: type3 has same/all nNA, same {min, q1}, smaller {max}, larger everything else, including sum
 # type3 vs type2: type3 has more/all nNA, smaller {q1, med, max}, larger {mean, q3, sum} ... much larger sum
 # combine_EDGAR_and_EPIC_emissions.r: combined (type3) output emissions=
+
+# type3 ~> EDGAR (but not much): from above:
+# combine_EDGAR_and_EPIC_emissions.r: EDGAR input emissions=
 #         cells=137241
 #         obs=137241
 #         min=0
 #         max=0.0861
 #         sum=264
 
-# type3 ~> EDGAR (but not much): from above:
-# combine_EDGAR_and_EPIC_emissions.r: EDGAR input emissions=
-#         cells=137241
-#         obs=137241
-#         min=0
-#         q1=0
-#         med=0.000157
-#         mean=0.00189
-#         q3=0.00312
-#         max=0.0861
-#         sum=259
-
 ### type4: use values from EPIC only where EPIC > EDGAR
 
-output.both.arr.type4 <- input.edgar.arr
-# we want to copy if either of
-# * is.na(EDGAR) & !is.na(EPIC)
-# * !is.na(EDGAR) & !is.na(EPIC) & (EPIC > EDGAR)
-# (EDGAR is never NA, but I'll test anyway, for safety)
-copy.criterion.type4 <-
-  (is.na(input.edgar.arr) & !is.na(input.epic.arr)) |
-  (!is.na(input.edgar.arr) & !is.na(input.epic.arr) & (input.epic.arr > input.edgar.arr))
-output.both.arr.type4[copy.criterion.type4] <- input.epic.arr[copy.criterion.type4]
+output.both.dv.arr.type4 <- input.epic.dv.arr # can't use EDGAR: dims don't match
+for (i.tstep in 1:output.both.dim.tsteps.n) {
+  for (i.lay in 1:output.both.dim.layers.n) {
+    input.epic.grid <- input.epic.dv.arr[,,i.lay,i.tstep]
+    output.grid <- input.edgar.arr
+    # we want to copy if either of
+    # * is.na(EDGAR) & !is.na(EPIC)
+    # * !is.na(EDGAR) & !is.na(EPIC) & (EPIC > EDGAR)
+    # (EDGAR-4.2 is never NA, but I'll test anyway, for safety)
+    copy.criterion.type4 <-
+      (is.na(input.edgar.arr) & !is.na(input.epic.grid)) |
+      (!is.na(input.edgar.arr) & !is.na(input.epic.grid) & (input.epic.grid > input.edgar.arr))
+    output.grid[copy.criterion.type4] <- input.epic.grid[copy.criterion.type4]
+    output.both.dv.arr.type4[,,i.lay,i.tstep] <- output.grid
+  }
+}
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(output.both.arr.type4)==', this_fn))
-print(dim(output.both.arr.type4))
+base::cat(base::sprintf('\n%s: dim(output.both.dv.arr.type4)==', this_fn))
+print(dim(output.both.dv.arr.type4))
 stats.to.stdout(
-  data=output.both.arr.type4,
+  data=output.both.dv.arr.type4,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: combined (type4) output emissions=', this_fn)
 )
 #         med=0.000297
 #         mean=0.00208
 #         q3=0.00341
-#         max=0.283
-#         sum=286
+#         max=0.0861
+#         sum=285
 
-# type4 ~> type3 ~> EDGAR (but not much more): from above:
+# type4 > type3 ~> EDGAR: from above:
 # combine_EDGAR_and_EPIC_emissions.r: EDGAR input emissions=
 #         cells=137241
 #         obs=137241
 
 #   end debugging-----------------------------------------------------
 
-### So I'll go with type4 for now:
-# output.both.arr <- output.both.arr.type4
-##### TODO: data is mirrored about horizontal axis! #####
-### so for now: reverse the _cols_ horizontally (not transpose)
-output.both.arr <-
-  output.both.arr.type4[,ncol(output.both.arr.type4):1] # flip the column indices
+### go with type4 for now:
+
+output.both.dv.arr <- output.both.dv.arr.type4
+##### TODO: determine why: data is mirrored about horizontal axis! #####
+# so for now: reverse the rows horizontally
+output.both.dv.arr <-   # reverse (only) the row indices
+  output.both.dv.arr.type4[,output.both.dim.y.n:1,,]
+
+# ----------------------------------------------------------------------
+# copy/mod IOAPI/fake datavar=TFLAG
+# ----------------------------------------------------------------------
+
+input.epic.TFLAG <- input.epic.fh$var[[input_epic_TFLAG_name]] # note double brackets!
+# input.epic.TFLAG.type <- input.epic.TFLAG$prec
+# arrggghhhh
+input.epic.TFLAG.type <- get.ncdf4.prec(input.epic.TFLAG)
+input.epic.TFLAG.mv <- input.epic.TFLAG$missval
+# get these from driver? no
+input.epic.TFLAG.units <- input.epic.TFLAG$units
+input.epic.TFLAG.long_name <- input.epic.TFLAG$longname
+# input.epic.TFLAG.var_desc <- input.epic.TFLAG$var_desc
+input.epic.TFLAG.var_desc <- ncdf4::ncatt_get(input.epic.fh, input_epic_TFLAG_name)$var_desc
+## and its dimensions? no, do that @ writetime, since dimensions are really about files
+
+output.both.TFLAG.var_desc <- input.epic.TFLAG.var_desc
 
 # ----------------------------------------------------------------------
 # define and write output
 # 6. (extra credit)remove output filehandle from the workspace/environment: call `rm(..)`
 
 ### 1. define output dimensions
+cat(sprintf('\n%s: 1. define output dimensions\n', this_fn)) # debugging
+
+## instead, copy input dimensions? I can't see how to :-(
+
+## 1.1. get values for output dimensions
+cat(sprintf('\n%s: 1.1. get values for output dimensions\n', this_fn)) # debugging
+
+# For non-grid dimensions, can use ncdf4
+
+input.epic.dim.tsteps.vals <- input.epic.dim.tsteps$vals
+# input.epic.dim.tsteps.units <- input.epic.dim.tsteps$units
+# ncdf4 (formatting added):
+# > [if] create_dimvar was specified to be FALSE, [then]
+# > the unit string MUST be empty ('') and
+# > the dimension values MUST be simple integers from 1 to the length of the dimension (e.g., 1:len)
+input.epic.dim.tsteps.units <- ''
+
+input.epic.dim.datetimes.vals <- input.epic.dim.datetimes$vals
+input.epic.dim.datetimes.units <- ''
+
+input.epic.dim.layers.vals <- input.epic.dim.layers$vals
+input.epic.dim.layers.units <- ''
+
+input.epic.dim.vars.vals <- input.epic.dim.vars$vals
+input.epic.dim.vars.units <- ''
+
+# For grid dimensions, I could use ncdf4, but, to overlay plots with "the usual" map,
+# we need dimensions uncentered and with units=m
+library(M3)
+input.epic.dim.x.vals <- M3::get.coord.for.dimension(
+  file=input_epic_fp, dimension="col", position="upper", units="m")$coords
+input.epic.dim.y.vals <- M3::get.coord.for.dimension(
+  file=input_epic_fp, dimension="row", position="upper", units="m")$coords
+
+## 1.2. Create output dimension definitions
+cat(sprintf('\n%s: 1.2. Create output dimension definitions\n', this_fn)) # debugging
+
+output.both.dim.x <- ncdim_def(
+  name=output_both_dim_x_name,
+  units=output_both_dim_x_units,
+  vals=input.epic.dim.x.vals,
+  unlim=FALSE,
+  create_dimvar=TRUE,
+# IOAPI does not create coordvars, but we will, for the grid (unless VERDI chokes)
+#   create_dimvar=FALSE,
+  calendar=NA,
+  longname=output_both_dim_x_long_name)
+
+output.both.dim.y <- ncdim_def(
+  name=output_both_dim_y_name,
+  units=output_both_dim_y_units,
+  vals=input.epic.dim.y.vals,
+  unlim=FALSE,
+  create_dimvar=TRUE,
+# IOAPI does not create coordvars, but we will, for the grid (unless VERDI chokes)
+#   create_dimvar=FALSE,
+  calendar=NA,
+  longname=output_both_dim_y_long_name)
+
+# output=input
+output.both.dim.datetimes <- ncdim_def(
+  name=input_epic_dim_datetimes_name,
+  units=input.epic.dim.datetimes.units,
+  vals=input.epic.dim.datetimes.vals,
+  unlim=FALSE,
+  create_dimvar=FALSE,
+  calendar=NA) # long_name omitted
+
+# output(LAY) == input(LAY)
+output.both.dim.layers <- ncdim_def(
+  name=input_epic_dim_layers_name,
+  units=input.epic.dim.layers.units,
+  vals=input.epic.dim.layers.vals,
+#  vals=c(1),
+#  vals=1,
+# dunno how this differs from the above, but below makes ncdf4 happy, and above did not
+#  vals=input.epic.dim.layers.vals[1],
+  unlim=FALSE,
+  create_dimvar=FALSE,
+  calendar=NA) # long_name omitted
+
+# output=input
+output.both.dim.tsteps <- ncdim_def(
+  name=input_epic_dim_tsteps_name,
+  units=input.epic.dim.tsteps.units,
+  vals=input.epic.dim.tsteps.vals,
+  unlim=FALSE,
+  create_dimvar=FALSE,
+  calendar=NA) # long_name omitted
+
+# output=input
+output.both.dim.vars <- ncdim_def(
+  name=input_epic_dim_vars_name,
+  units=input.epic.dim.vars.units,
+  vals=input.epic.dim.vars.vals,
+  unlim=FALSE,
+  create_dimvar=FALSE,
+  calendar=NA) # long_name omitted
 
-## instead, copy input dimensions? No: above have type=list, not type=ncdim4 :-(
+### 2. define output datavar(s)
+cat(sprintf('\n%s: 2. define output datavar(s)\n', this_fn)) # debugging
+
+## Try to make TFLAG the first "datavar"? Not sure it matters, but in the IOAPI files I've `ncdump`ed, it always is.
+
+output.both.TFLAG <- ncdf4::ncvar_def(
+  name=output_both_TFLAG_name,
+  units=input.epic.TFLAG.units,
+# (TSTEP, VAR, DATE-TIME) is how it looks in `ncdump`, but ncdf4 wants ...
+  dim=list(output.both.dim.datetimes, output.both.dim.vars, output.both.dim.tsteps),
+  prec=input.epic.TFLAG.type,
+# ncdf4::ncatt_put vomits on being asked to create datavar of type=integer with a missing_value:
+# > [1] "ncatt_put_inner: prec to create: integer"
+# > Error in R_nc4_put_att_logical: asked to put a NA value, but the variable's type is not a double or float, which are the only two types that have a defined NaN value
+# > [1] "Error in ncatt_put, while writing attribute _FillValue with value NA"
+# > Error in ncatt_put_inner(ncid2use, newvar$id, "_FillValue", v$missval,  : 
+# >   Error return from C call R_nc4_put_att_logical for attribute _FillValue
+# > Calls: <Anonymous> -> ncvar_add -> ncatt_put_inner
+# so don't try!
+#   missval=input.epic.TFLAG.mv,
+  longname=input.epic.TFLAG.long_name,
+  shuffle=FALSE,
+  compression=NA,
+  chunksizes=NA)
+
+output.both.dv <- ncdf4::ncvar_def(
+  name=output_both_dv_name,
+  units=output_both_dv_units,
+# (TSTEP, LAY, ROW, COL) is how it looks in `ncdump`, but ncdf4 wants ...
+  dim=list(output.both.dim.x, output.both.dim.y, output.both.dim.layers, output.both.dim.tsteps),
+  missval=input.epic.dv.mv,
+  longname=output_both_dv_long_name,
+  prec=input.epic.dv.type,
+  shuffle=FALSE,
+  compression=NA,
+  chunksizes=NA)
 
 input.epic.dim.x.vals <- input.epic.dim.x$vals
 input.epic.dim.y.vals <- input.epic.dim.y$vals
   longname=output_both_dim_y_long_name)
 
 ### 2. define output datavar(s)
+cat(sprintf('\n%s: 2. define output datavar(s)\n', this_fn)) # debugging
+
+## Try to make TFLAG the first "datavar"? Not sure it matters, but in the IOAPI files I've `ncdump`ed, it always is.
+
+output.both.TFLAG <- ncdf4::ncvar_def(
+  name=output_both_TFLAG_name,
+  units=input.epic.TFLAG.units,
+# (TSTEP, VAR, DATE-TIME) is how it looks in `ncdump`, but ncdf4 wants ...
+  dim=list(output.both.dim.datetimes, output.both.dim.vars, output.both.dim.tsteps),
+  prec=input.epic.TFLAG.type,
+# ncdf4::ncatt_put vomits on being asked to create datavar of type=integer with a missing_value:
+# > [1] "ncatt_put_inner: prec to create: integer"
+# > Error in R_nc4_put_att_logical: asked to put a NA value, but the variable's type is not a double or float, which are the only two types that have a defined NaN value
+# > [1] "Error in ncatt_put, while writing attribute _FillValue with value NA"
+# > Error in ncatt_put_inner(ncid2use, newvar$id, "_FillValue", v$missval,  : 
+# >   Error return from C call R_nc4_put_att_logical for attribute _FillValue
+# > Calls: <Anonymous> -> ncvar_add -> ncatt_put_inner
+# so don't try!
+#   missval=input.epic.TFLAG.mv,
+  longname=input.epic.TFLAG.long_name,
+  shuffle=FALSE,
+  compression=NA,
+  chunksizes=NA)
 
 output.both.dv <- ncdf4::ncvar_def(
   name=output_both_dv_name,
   units=output_both_dv_units,
 #  dim=list(output.both.dim.y, output.both.dim.x),
-# (ROW, COL) is how it looks in `ncdump`, but internally we want ...
-  dim=list(output.both.dim.x, output.both.dim.y),
+# (TSTEP, LAY, ROW, COL) is how it looks in `ncdump`, but ncdf4 wants the reverse
+  dim=list(output.both.dim.x, output.both.dim.y, output.both.dim.layers, output.both.dim.tsteps),
   missval=input.epic.dv.mv,
   longname=output_both_dv_long_name,
   prec=input.epic.dv.type,
   chunksizes=NA)
 
 ### 3. create output file
+cat(sprintf('\n%s: 3. create output file\n', this_fn)) # debugging
 
 # # start debugging-----------------------------------------------------
 # base::cat(base::sprintf('\n%s: about to attempt\n', this_fn))
 
 output.both.fh <- ncdf4::nc_create(
   filename=output_both_fp,
-  vars=output.both.dv,               # can also be vector or list
+#   vars=output.both.dv,               # can also be vector or list? no, just a list
+# Try to make TFLAG the first "datavar"? Not sure it matters, but in the IOAPI files I've `ncdump`ed, it always is.
+  vars=list(output.both.TFLAG, output.both.dv),
   force_v4=FALSE,
 # # to debug, toggle verbose
 #   verbose=TRUE)
   verbose=FALSE)
 
 ## 3.1. create output file/global attributes
+cat(sprintf('\n%s: 3.1. create output file/global attributes\n', this_fn)) # debugging
 
-## copy all attributes from EPIC input
+## copy all attributes from EPIC input, except :FILEDESC
 ## (see extended discussion in comments to compute_EPIC_emissions.r)
 
-input.epic.attrs <- ncdf4::ncatt_get(nc=input.epic.fh, varid=0,
-  attname=NA) # this is what gets us a list of all file/global attributes
+input.epic.attrs <- ncdf4::ncatt_get(nc=input.epic.fh, 
+  varid=0, attname=NA) # this is what gets us a list of all file/global attributes
 input.epic.attrs.n <- length(input.epic.attrs)
 if (input.epic.attrs.n > 0) {
   input.epic.attrs.names <- names(input.epic.attrs)
   for (i.attrs in 1:input.epic.attrs.n) {
     attr.name <- input.epic.attrs.names[i.attrs]
-    attr.val <- input.epic.attrs[[i.attrs]] # copy from input
-    attr.type <- typeof(attr.val)
+    if        (attr.name == 'FILEDESC') {
+      attr.val <- output_both_attr_filedesc
+      attr.type <- output_both_attr_filedesc_type
+    } else {
+      attr.val <- input.epic.attrs[[i.attrs]] # copy from input
+      attr.type <- typeof(attr.val)
+    }
     ncdf4::ncatt_put(
       nc=output.both.fh,
       varid=0,            # for file/global attributes
 } # end if (input.epic.attrs.n > 0)
 
 ### 4. write data to output datavar
+cat(sprintf('\n%s: 4. write data to output datavar\n', this_fn)) # debugging
 
 ## 4.1. write datavar data
+cat(sprintf('\n%s: 4.1. write datavar data\n', this_fn)) # debugging
 
 # # start debugging-----------------------------------------------------
 # base::cat(base::sprintf('\n%s: about to attempt\n', this_fn))
 # base::cat(base::sprintf('ncdf4::ncvar_put(\n'))
 # base::cat(base::sprintf('  nc=output.both.fh == %s\n', output_both_fp))
 # base::cat(base::sprintf('  varid=%s\n', output_both_dv_name))
-# base::cat(base::sprintf('  vals=output.both.arr == '))
-# print(dim(output.both.arr))
+# base::cat(base::sprintf('  vals=output.both.dv.arr == '))
+# print(dim(output.both.dv.arr))
 # #   end debugging-----------------------------------------------------
 
 ncdf4::ncvar_put(
   nc=output.both.fh,
+#   varid=output_both_TFLAG_name,
+  varid=output.both.TFLAG,
+  vals=output.both.TFLAG.arr,
+  start=output.both.TFLAG.write.start,
+  count=output.both.TFLAG.write.count,
+# # to debug, toggle verbose
+#   verbose=TRUE)
+  verbose=FALSE)
+
+# ======================================================================
+
+ncdf4::ncvar_put(
+  nc=output.both.fh,
 #  varid=output.both.dv,
   varid=output_both_dv_name,
-  vals=output.both.arr,
-  start=output.both.write.start,
-  count=output.both.write.count,
+  vals=output.both.dv.arr,
+  start=output.both.dv.write.start,
+  count=output.both.dv.write.count,
 # # to debug, toggle verbose
 #   verbose=TRUE)
   verbose=FALSE)
 
-## 4.1. write datavar metadata (i.e., additional attributes)
+## 4.2. write datavar metadata (i.e., additional attributes)
+cat(sprintf('\n%s: 4.2. write datavar metadata (i.e., additional attributes)\n', this_fn)) # debugging
+
+ncdf4::ncatt_put(
+  nc=output.both.fh,
+  varid=output.both.TFLAG,
+  attname='var_desc',
+  attval=output.both.TFLAG.var_desc,
+  prec='text',
+  definemode=FALSE,
+# # to debug, toggle verbose
+#   verbose=TRUE)
+  verbose=FALSE)
 
 ncdf4::ncatt_put(
   nc=output.both.fh,
 rm(input.edgar.fh)
 rm(input.epic.fh)
 
+#----------------------------------------------------------------------
+# verify VERDI compatibility? output=combined emissions
+#----------------------------------------------------------------------
+
+  ### TODO: integrate VERDI into bash_utilities.sh
+  ### plot our datavar of interest with VERDI
+  output_both_dv_name_VERDI <- sprintf('%s[1]', output_both_dv_name)
+  ## VERDI really wants an FQP
+  output_both_fp_VERDI <- system(sprintf('readlink -f %s', output_both_fp), intern=TRUE)
+  output_both_dv_VERDI_cmd <-
+    sprintf('verdi -f %s -s %s -gtype tile &',
+      output_both_fp_VERDI, output_both_dv_name_VERDI)
+  cat(sprintf("\n%s: debugging IOAPI: about to run '%s'\n", this_fn, output_both_dv_VERDI_cmd))
+#   cat(sprintf("\n%s: not debugging IOAPI: skipping '%s'\n", this_fn, output_both_dv_VERDI_cmd))
+  system(output_both_dv_VERDI_cmd)
+
 # ----------------------------------------------------------------------
 # visualize output emissions
 # ----------------------------------------------------------------------
 # something whacks the CRS, so
 output.both.raster@crs <- aqme.crs
 
-# start debugging-----------------------------------------------------
-cat('\noutput.both.raster==\n')
-output.both.raster
-# class       : RasterLayer 
-# dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
-# resolution  : 12000, 12000  (x, y)
-# extent      : -2550000, 2958000, -1722000, 1866000  (xmin, xmax, ymin, ymax)
-# coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 
-# data source : .../5yravg_20111219_pure_emis_both.nc 
-# names       : N2O 
-# zvar        : N2O 
-
-cat('output.both.raster@extent==\n')
-output.both.raster@extent
-# class       : Extent 
-# xmin        : -2550000 
-# xmax        : 2958000 
-# ymin        : -1722000 
-# ymax        : 1866000 
-#   end debugging-----------------------------------------------------
+# # start debugging-----------------------------------------------------
+# cat('\noutput.both.raster==\n')
+# output.both.raster
+# # class       : RasterLayer 
+# # dimensions  : 299, 459, 137241  (nrow, ncol, ncell)
+# # resolution  : 12000, 12000  (x, y)
+# # extent      : -2550000, 2958000, -1722000, 1866000  (xmin, xmax, ymin, ymax)
+# # coord. ref. : +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000 
+# # data source : .../5yravg_20111219_pure_emis_both.nc 
+# # names       : N2O 
+# # zvar        : N2O 
+
+# cat('output.both.raster@extent==\n')
+# output.both.raster@extent
+# # class       : Extent 
+# # xmin        : -2550000 
+# # xmax        : 2958000 
+# # ymin        : -1722000 
+# # ymax        : 1866000 
+# #   end debugging-----------------------------------------------------
 
 # TODO: if using output as args below, `project.NorAm.boundaries.for.CMAQ` fails with
 # > Grid type 0 cannot be handled by this function.

File compute_EPIC_emissions.r

 ## EPIC N2O emissions over AQMEII-NA (actually CONUS): reunit-ed, monotonic, VERDI-compliant
 
 output_epic_fp <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_FP', unset=NA) # path to output file
+
 # file/global attributes
-output_epic_attr_nlays = as.numeric(Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS'))
-output_epic_attr_nlays_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS_TYPE')
-output_epic_attr_varlist = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST')
-output_epic_attr_varlist_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST_TYPE')
-output_epic_attr_filedesc = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC')
-output_epic_attr_filedesc_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC_TYPE')
+output_epic_attr_nlays = as.numeric(Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS', unset=NA))
+output_epic_attr_nlays_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_NLAYS_TYPE', unset=NA)
+output_epic_attr_varlist = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST', unset=NA)
+output_epic_attr_varlist_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_VARLIST_TYPE', unset=NA)
+output_epic_attr_filedesc = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC', unset=NA)
+output_epic_attr_filedesc_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC_TYPE', unset=NA)
+
 # dimensions
+output_epic_dim_datetimes_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_NAME', unset=NA)
+output_epic_dim_tsteps_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_NAME', unset=NA)
+output_epic_dim_vars_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_VAR_NAME', unset=NA)
 output_epic_dim_x_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_NAME', unset=NA)
 output_epic_dim_x_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_UNITS', unset=NA)
 output_epic_dim_x_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_LONG_NAME', unset=NA)
 output_epic_dim_layers_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_UNITS', unset=NA)
 output_epic_dim_layers_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_LONG_NAME', unset=NA)
 
-output_epic_dim_datetimes_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_NAME', unset=NA)
-output_epic_dim_tsteps_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_NAME', unset=NA)
-output_epic_dim_vars_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_VAR_NAME', unset=NA)
-
 # datavars
 output_epic_dv_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_NAME', unset=NA)
 output_epic_dv_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_UNITS', unset=NA)
 # internal functions (defined)
 # ----------------------------------------------------------------------
 
+### TODO: move to regrid_utils: get.ncdf4.prec also used by combine_EDGAR_and_EPIC_emissions.r
+
 # workaround apparent bug in ncdf4
 get.ncdf4.prec <- function(
   ncvar4 # object of type=ncvar4
 
 ### "fake" datavar=TFLAG, for IOAPI and (IIUC) VERDI
 input.epic.TFLAG.arr <- ncdf4::ncvar_get(input.epic.fh, varid=input_epic_TFLAG_name)
-input.epic.TFLAG.dims <- dim(input.epic.TFLAG.arr)
-input.epic.TFLAG.dims.n <- length(input.epic.TFLAG.dims)
 input.epic.TFLAG <- input.epic.fh$var[[input_epic_TFLAG_name]]
 input.epic.TFLAG.type <- input.epic.TFLAG$prec
 input.epic.TFLAG.mv <- input.epic.TFLAG$missval
 
-### get input dimensions useful for output (i.e., of type=ncdf4, not type=list)
-# # print(input.epic.fh$dim)
-# output.epic.dim.y <- input.epic.fh$dim[2]
-# output.epic.dim.x <- input.epic.fh$dim[3]
-input.epic.dim.y <- input.epic.fh$dim[[input_epic_dim_y_name]]
-input.epic.dim.x <- input.epic.fh$dim[[input_epic_dim_x_name]]
-input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
-input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
-input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
-input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
-
 # ----------------------------------------------------------------------
 # get input and output dimensionalities
 # ----------------------------------------------------------------------
 ### Both EPIC input and output have same spatial grid, input also has layers=crops as leftmost dim (reversed by R)
 # input.epic.dv.dims <- dim(input.epic.dv)
 # > NULL
-input.epic.dv.dims <- dim(input.epic.dv.arr)
-input.epic.dv.dims.n <- length(input.epic.dv.dims)
 
-input.epic.dim.x.n <- input.epic.dim.x$len
-output.epic.dim.x.n <- input.epic.dim.x.n
+### EPIC input (emissivities)
 
-input.epic.dim.y.n <- input.epic.dim.y$len
-output.epic.dim.y.n <- input.epic.dim.y.n
+## EPIC input dims
 
-input.epic.dim.layers.n <- input.epic.dim.layers$len
-output.epic.dim.layers.n <- 1 # in output, `layer`=verticality, !=crop
+input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
+input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
+input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
+input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
+input.epic.dim.x <- input.epic.fh$dim[[input_epic_dim_x_name]]
+input.epic.dim.y <- input.epic.fh$dim[[input_epic_dim_y_name]]
 
+input.epic.dim.datetimes.n <- input.epic.dim.datetimes$len
+input.epic.dim.layers.n <- input.epic.dim.layers$len
 input.epic.dim.tsteps.n <- input.epic.dim.tsteps$len
-output.epic.dim.tsteps.n <- input.epic.dim.tsteps.n
+input.epic.dim.vars.n <- input.epic.dim.vars$len
+input.epic.dim.x.n <- input.epic.dim.x$len
+input.epic.dim.y.n <- input.epic.dim.y$len
 
+## EPIC input "real" datavar
+
+input.epic.dv.dims <- dim(input.epic.dv.arr)
+input.epic.dv.dims.n <- length(input.epic.dv.dims)
 # spatial cells only: ignore layers, tsteps
 input.epic.dv.cells.n <- input.epic.dim.y.n * input.epic.dim.x.n
-output.epic.dv.cells.n <- input.epic.dv.cells.n
 
-input.epic.dim.datetimes.n <- input.epic.dim.datetimes$len
-output.epic.dim.datetimes.n <- input.epic.dim.datetimes.n
+## EPIC input TFLAG
 
-input.epic.dim.vars.n <- input.epic.dim.vars$len
+input.epic.TFLAG.dims <- dim(input.epic.TFLAG.arr)
+input.epic.TFLAG.dims.n <- length(input.epic.TFLAG.dims)
+
+### EPIC output (emissions)
+
+output.epic.dim.datetimes.n <- input.epic.dim.datetimes.n
+output.epic.dim.layers.n <- 1 # in output, `layer`=verticality, !=crop
+output.epic.dim.tsteps.n <- input.epic.dim.tsteps.n
 output.epic.dim.vars.n <- input.epic.dim.vars.n
+output.epic.dim.x.n <- input.epic.dim.x.n
+output.epic.dim.y.n <- input.epic.dim.y.n
 
-# output.epic.dv.dims <- c(output.epic.dim.y.n, output.epic.dim.x.n)
-# (row, col) is how we see it, but
-# > dim(input.epic.dv.arr)
-# [1] 459 299  42
-# output.epic.dv.dims <- c(output.epic.dim.x.n, output.epic.dim.y.n)
-# Also, gotta be VERDI-compliant
+output.epic.dv.cells.n <- input.epic.dv.cells.n
 output.epic.dv.dims <-
   c(output.epic.dim.x.n, output.epic.dim.y.n, output.epic.dim.layers.n, output.epic.dim.tsteps.n)
 output.epic.dv.dims.n <- length(output.epic.dv.dims)
 ## other than the implicit/degenerate model year ... which, for VERDI-compliance, it appears we must use :-(
 ## So gotta deal with implicit/degenerate TSTEP==1 ...
 
-### Still need to set read start point and read count (both by datavar dimensions)?
 input.epic.dv.read.start <- rep(1, input.epic.dv.dims.n)   # start=(1,1,1,...)
 input.epic.dv.read.count <- input.epic.dv.dims             # count=(n,n,n,...)
 
 
 } else {
   # TODO: throw
-#   cat(sprintf('\n%s: ERROR: dimensionality=%i < 3 of EPIC input datavar=%s\n',
-#     this_fn, input.epic.dv.dims.n, input_epic_dv_name))
   stop(sprintf('%s: ERROR: dimensionality=%i > 4 of EPIC input datavar=%s',
     this_fn, input.epic.dv.dims.n, input_epic_dv_name))
 } # end testing input.epic.dv.dims.n
 #   sig.digs=sigdigs,
 #   title=base::sprintf('\n%s: EPIC input emittivities=', this_fn)
 # )
+
 # # compute_EPIC_emissions.r: EPIC input emittivities=
 # #         cells=5764122
 # #         obs=369846
 
 # For non-grid dimensions, can use ncdf4
 
-input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
+# input.epic.dim.tsteps <- input.epic.fh$dim[[input_epic_dim_tsteps_name]]
 input.epic.dim.tsteps.vals <- input.epic.dim.tsteps$vals
 # input.epic.dim.tsteps.units <- input.epic.dim.tsteps$units
 # ncdf4 (formatting added):
 # > the dimension values MUST be simple integers from 1 to the length of the dimension (e.g., 1:len)
 input.epic.dim.tsteps.units <- ''
 
-input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
+# input.epic.dim.datetimes <- input.epic.fh$dim[[input_epic_dim_datetimes_name]]
 input.epic.dim.datetimes.vals <- input.epic.dim.datetimes$vals
 # input.epic.dim.datetimes.units <- input.epic.dim.datetimes$units
 input.epic.dim.datetimes.units <- ''
 
-input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
+# input.epic.dim.layers <- input.epic.fh$dim[[input_epic_dim_layers_name]]
 input.epic.dim.layers.vals <- input.epic.dim.layers$vals
 # input.epic.dim.layers.units <- input.epic.dim.layers$units
 input.epic.dim.layers.units <- ''
 
-## input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
+# input.epic.dim.vars <- input.epic.fh$dim[[input_epic_dim_vars_name]]
 input.epic.dim.vars.vals <- input.epic.dim.vars$vals
 # input.epic.dim.vars.units <- input.epic.dim.vars$units
 input.epic.dim.vars.units <- ''