Commits

Tom Roche committed 6d55d0c

(INTERIM COMMIT, BREAKS LATER CODE) steps 1-2 of VERDIfication

EJC (et probably al) needs VERDI for visualization.
VERDI requires IOAPI-compliance (to a high degree, anyway).
So VERDIfied the following (e.g., adding "datavar"=TFLAG, adding dims={DATE-TIME, LAY, TSTEP, VAR),
with VERDIfication verified by actually running VERDI from the script(s) against the following file/datavar tuples on terrae:

* EPIC raw emittivities, DN2: just verified
* EPIC pure emittivities, N2Oemittivities: VERDIfied
* EPIC emissions, N2O: VERDIfied

Did not do the EDGAR emissions, have not (but must) do the combined emissions!
TODO: combine_EDGAR_and_EPIC_emissions.r: currently broken by these changes!
Also:

* 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), but combine_EDGAR_and_EPIC_emissions.r is currently broken by these changes!

TODO:
* handle {name, path to} VERDI in regrid_utils/bash_utilities.sh
*** make combine(EDGAR, EPIC) VERDI-loadable!
* 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)

Comments (0)

Files changed (4)

 
 ## datavar dims and attributes
 # TODO: get from template file or from EPIC
+# "real" data
 AQMEIINA_DV_NAME='N2O'
 # IOAPI pads varattr=long_name to length=16 with trailing spaces
 AQMEIINA_DV_LONG_NAME="$(printf '%-16s' ${AQMEIINA_DV_NAME})"
 # Don't single-quote the payload: double-quote it (OK inside parens inside double-quotes)
 AQMEIINA_DV_VAR_DESC="$(printf '%-80s' "Model species ${AQMEIINA_DV_NAME}")"
 
+# "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
+AQMEIINA_TFLAG_NAME='TFLAG'
+# 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>')"
+# 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_DIM_LAYER_NAME='LAY'
 AQMEIINA_DIM_LAYER_LONG_NAME='index of layers above surface'
 AQMEIINA_DIM_LAYER_UNITS='unitless'
 
-AQMEIINA_DIM_TIME_NAME='TSTEP'
-# AQMEIINA_DIM_TIME_UNITS # vary
-# AQMEIINA_DIM_TIME_LONG_NAME='timestep'
+AQMEIINA_DIM_TSTEP_NAME='TSTEP'
+AQMEIINA_DIM_TSTEP_UNITS='' # they vary!
+AQMEIINA_DIM_TSTEP_LONG_NAME='timestep'
 
 AQMEIINA_DIM_X_N='459'
 AQMEIINA_DIM_X_NAME='COL'
 # AQMEIINA_DIM_Y_LONG_NAME='Fortran-style index to grid columns from lower-left origin' # my invention (TODO: CHECK it's not from top)
 AQMEIINA_DIM_Y_LONG_NAME="${AQMEIINA_DIM_X_LONG_NAME}"
 
+# dimensions we don't really need, but VERDI does
+AQMEIINA_DIM_DATETIME_NAME='DATE-TIME'
+AQMEIINA_DIM_LAYER_NAME='LAY'
+AQMEIINA_DIM_TSTEP_NAME='TSTEP'
+AQMEIINA_DIM_VAR_NAME='VAR'
+
 ### artifact visualization
 
 ## for visualization (generally)
 export AQMEII_EDGAR_N2O_REGRID_DIM_LAYER_LONG_NAME="${AQMEIINA_DIM_LAYER_LONG_NAME}"
 export AQMEII_EDGAR_N2O_REGRID_DIM_LAYER_N='1'
 
-export AQMEII_EDGAR_N2O_REGRID_DIM_TIME_NAME="${AQMEIINA_DIM_TIME_NAME}"
-export AQMEII_EDGAR_N2O_REGRID_DIM_TIME_UNITS="${AQMEIINA_DIM_TIME_UNITS}"
-export AQMEII_EDGAR_N2O_REGRID_DIM_TIME_LONG_NAME="${AQMEIINA_DIM_TIME_LONG_NAME}"
-export AQMEII_EDGAR_N2O_REGRID_DIM_TIME_N='1'
+export AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_NAME="${AQMEIINA_DIM_TSTEP_NAME}"
+export AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_UNITS="${AQMEIINA_DIM_TSTEP_UNITS}"
+export AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_LONG_NAME="${AQMEIINA_DIM_TSTEP_LONG_NAME}"
+export AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_N='1'
 
 export AQMEII_EDGAR_N2O_REGRID_DIM_X_NAME="${AQMEIINA_DIM_X_NAME}"
 export AQMEII_EDGAR_N2O_REGRID_DIM_X_UNITS="${AQMEIINA_DIM_X_UNITS}"
 export AQMEII_EPIC_N2O_EMITT_DIM_X_NAME="${AQMEIINA_DIM_X_NAME}"
 export AQMEII_EPIC_N2O_EMITT_DIM_Y_NAME="${AQMEIINA_DIM_Y_NAME}"
 export AQMEII_EPIC_N2O_EMITT_DIM_LAYERS_NAME="${AQMEIINA_DIM_LAYER_NAME}"
+export AQMEII_EPIC_N2O_EMITT_DIM_TSTEPS_NAME="${AQMEIINA_DIM_TSTEP_NAME}"
 
-## datavar
+## "real" datavar, with information we actually want
 export AQMEII_EPIC_N2O_EMITT_DATAVAR_NAME='N2Oemittivities' # dim=(LAY, ROW, COL)
 # IOAPI pads varattr=long_name to length=16 with trailing spaces
 export AQMEII_EPIC_N2O_EMITT_DATAVAR_LONG_NAME="$(printf '%-16s' 'N2O emittivities')"
 # IOAPI pads varattr=var_desc to length=80 with trailing spaces
 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'
+# 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'
+
+# get these from driver
+# input.epic.TFLAG.units <- input.epic.TFLAG$units
+# input.epic.TFLAG.long_name <- input.epic.TFLAG$long_name
+# input.epic.TFLAG.var_desc <- input.epic.TFLAG$var_desc
+# input_epic_dim_datetimes_name]]
+# input_epic_dim_layers_name]]
+# input_epic_dim_vars_name]]
+
 ### intermediate product: EPIC N2O emissions over AQMEII in mol/s
 
 AQMEII_EPIC_N2O_EMISS_ROOT="${AQMEII_EPIC_N2O_EMITT_ROOT}_emis"
 ## dimensions
 # names: unfortunately ncdf4 appears to need these defined to walk its ADSs :-(
 export AQMEII_EPIC_N2O_EMISS_DIM_X_NAME="${AQMEIINA_DIM_X_NAME}"
-export AQMEII_EPIC_N2O_EMISS_DIM_Y_NAME="${AQMEIINA_DIM_Y_NAME}"
 # might as well make these explicit, if we're gonna create coordvars
 export AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_UNITS="${AQMEIINA_DIM_X_UNITS}"
 export AQMEII_EPIC_N2O_EMISS_DIM_X_ATTR_LONG_NAME="${AQMEIINA_DIM_X_LONG_NAME}"
+
+export AQMEII_EPIC_N2O_EMISS_DIM_Y_NAME="${AQMEIINA_DIM_Y_NAME}"
 export AQMEII_EPIC_N2O_EMISS_DIM_Y_ATTR_UNITS="${AQMEIINA_DIM_Y_UNITS}"
 export AQMEII_EPIC_N2O_EMISS_DIM_Y_ATTR_LONG_NAME="${AQMEIINA_DIM_Y_LONG_NAME}"
 
+export AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_NAME="${AQMEIINA_DIM_DATETIME_NAME}"
+export AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_ATTR_UNITS="${AQMEIINA_DIM_DATETIME_UNITS}"
+export AQMEII_EPIC_N2O_EMISS_DIM_DATETIME_ATTR_LONG_NAME="${AQMEIINA_DIM_DATETIME_LONG_NAME}"
+
+export AQMEII_EPIC_N2O_EMISS_DIM_LAYER_NAME="${AQMEIINA_DIM_LAYER_NAME}"
+export AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_UNITS="${AQMEIINA_DIM_LAYER_UNITS}"
+export AQMEII_EPIC_N2O_EMISS_DIM_LAYER_ATTR_LONG_NAME="${AQMEIINA_DIM_LAYER_LONG_NAME}"
+
+export AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_NAME="${AQMEIINA_DIM_TSTEP_NAME}"
+export AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_ATTR_UNITS="${AQMEIINA_DIM_TSTEP_UNITS}"
+export AQMEII_EPIC_N2O_EMISS_DIM_TSTEP_ATTR_LONG_NAME="${AQMEIINA_DIM_TSTEP_LONG_NAME}"
+
+export AQMEII_EPIC_N2O_EMISS_DIM_VAR_NAME="${AQMEIINA_DIM_VAR_NAME}"
+export AQMEII_EPIC_N2O_EMISS_DIM_VAR_ATTR_UNITS="${AQMEIINA_DIM_VAR_UNITS}"
+export AQMEII_EPIC_N2O_EMISS_DIM_VAR_ATTR_LONG_NAME="${AQMEIINA_DIM_VAR_LONG_NAME}"
+
 ## datavar
-export AQMEII_EPIC_N2O_EMISS_DATAVAR_NAME="${AQMEIINA_DV_NAME}" # dim=(ROW, COL)
+export AQMEII_EPIC_N2O_EMISS_DATAVAR_NAME="${AQMEIINA_DV_NAME}"
 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}"
+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 GFED_REGRID_HOURLY_EMIS_N_LAYER='1' # one species in file
 
 # # datavar dims and attributes
-# export GFED_REGRID_HOURLY_EMIS_DIM_TIME_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TIME_NAME}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_TIME_UNITS="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TIME_UNITS}"
-# export GFED_REGRID_HOURLY_EMIS_DIM_TIME_LONG_NAME="${GFED_GLOBAL_3HOURLY_FRAC_CONV_DIM_TIME_LONG_NAME}"
+# 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}"
   'BELD_CSV_to_RDS' \
   'distill_EPIC' \
   'compute_EPIC' \
-  'combine_emis' \
   'teardown' \
 ; do
   echo -e "\n$ ${THIS_FN}::main loop::${CMD}\n"

compute_EPIC_emissions.r

 
 beld_fp <- Sys.getenv('BELD_RDS_FP', unset=NA) # path to serialized BELD array
 
-## "distilled" EPIC N2O emittivities: reunit-ed, monotonic
+## "distilled" EPIC N2O emittivities: reunit-ed, monotonic, VERDI-compliant
 
 input_epic_fp <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_FP', unset=NA) # path to input file
 input_epic_dv_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_DATAVAR_NAME', unset=NA)
+input_epic_dim_tsteps_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMITT_DIM_TSTEPS_NAME', unset=NA)
 input_epic_dim_layers_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMITT_DIM_LAYERS_NAME', unset=NA)
 input_epic_dim_x_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_DIM_X_NAME', unset=NA)
 input_epic_dim_y_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_DIM_Y_NAME', unset=NA)
 
+# "fake" datavar=TFLAG, required by IOAPI (and more to the point (IIUC), VERDI)
+input_epic_TFLAG_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_TFLAG_NAME', unset=NA)
+input_epic_dim_datetimes_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_DATETIME_NAME', unset=NA)
+input_epic_dim_vars_name <- Sys.getenv('AQMEII_EPIC_N2O_EMITT_TFLAG_DIM_VAR_NAME', unset=NA)
+
 ## outputs
 
-## EPIC N2O emissions
+## 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_filedesc = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC')
 output_epic_attr_filedesc_type = Sys.getenv('AQMEII_EPIC_N2O_EMISS_ATTR_FILEDESC_TYPE')
 # dimensions
-output_epic_dim_layers_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_NAME', unset=NA)
-output_epic_dim_layers_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_ATTR_UNITS', unset=NA)
-output_epic_dim_layers_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYERS_ATTR_LONG_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_y_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_Y_NAME', unset=NA)
 output_epic_dim_y_units <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_Y_ATTR_UNITS', unset=NA)
 output_epic_dim_y_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_Y_ATTR_LONG_NAME', unset=NA)
+output_epic_dim_layers_name <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_DIM_LAYER_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)
 output_epic_dv_long_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_DATAVAR_ATTR_LONG_NAME', unset=NA)
 output_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)
+output_epic_TFLAG_name <- Sys.getenv('AQMEII_EPIC_N2O_EMISS_TFLAG_NAME', unset=NA)
+# output_epic_TFLAG_var_desc <-  Sys.getenv('AQMEII_EPIC_N2O_EMISS_TFLAG_ATTR_VAR_DESC', unset=NA) # get/copy programmatically
 
 # visualization
 
 # internal functions (defined)
 # ----------------------------------------------------------------------
 
+# 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)
+  }
+}
+
 ### Return TRUE if value is NA or 0, else FALSE. FOR USE ONLY IN `if` CONDITIONS (IIRC)!
 ### ASSERT: input is numeric (we won't check!)
 is.val.na.or.zero <- function(val) {
 
 library(ncdf4)
 
-### EPIC per-crop N2O emittivities, previously "distilled"
 input.epic.fh <- ncdf4::nc_open(input_epic_fp, write=FALSE, readunlim=FALSE)
+
+### EPIC per-crop N2O emittivities, previously "distilled"
 # 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.layers <- input.epic.fh$dim[[input_epic_dim_layers_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.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.arr)
+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
 input.epic.dim.layers.n <- input.epic.dim.layers$len
 output.epic.dim.layers.n <- 1 # in output, `layer`=verticality, !=crop
 
-input.epic.dv.cells.n <- input.epic.dim.y.n * input.epic.dim.x.n # spatial only
-output.epic.dv.cells.n <- input.epic.dv.cells.n # same spatiality
+input.epic.dim.tsteps.n <- input.epic.dim.tsteps$len
+output.epic.dim.tsteps.n <- input.epic.dim.tsteps.n
+
+# 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
+
+input.epic.dim.vars.n <- input.epic.dim.vars$len
+output.epic.dim.vars.n <- input.epic.dim.vars.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.arr)
+# > dim(input.epic.dv.arr)
 # [1] 459 299  42
-output.epic.dv.dims <- c(output.epic.dim.x.n, output.epic.dim.y.n)
+# output.epic.dv.dims <- c(output.epic.dim.x.n, output.epic.dim.y.n)
+# Also, gotta be VERDI-compliant
+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)
 
 # ----------------------------------------------------------------------
 
 ### EPIC input
 
-## TODO: "Pierce-style read" (see help(ncvar_get)#Examples): protects against data overload (somewhat) but more complicated to implement.
-## Normally one reads one timestep at a time, but EPIC input has no timesteps (other than the implicit/degenerate model year).
+## 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 TSTEP==1 ...
 
 ### Still need to set read start point and read count (both by datavar dimensions)?
-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.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,...)
+
+# so gotta handle TSTEP=1: remembering timelike dim is always the LAST dimension!
+if      (input.epic.dv.dims.n < 3) {
+  # 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 < 3 of EPIC input datavar=%s',
+    this_fn, input.epic.dv.dims.n, input_epic_dv_name))
+
+} else if (input.epic.dv.dims.n == 3) {
+  if (input.epic.dim.tsteps.n != 1) {
+    stop(sprintf('%s: ERROR: dimensionality=3 of EPIC input datavar=%s, but dim(tsteps)=%i',
+      this_fn, input_epic_dv_name, input.epic.dim.tsteps.n))
+  }
+  # append tsteps
+  input.epic.dv.dims <- c(input.epic.dv.dims, 1)
+  input.epic.dv.read.count <- input.epic.dv.dims
+  input.epic.dv.read.start <- c(input.epic.dv.read.start, 1)
+  input.epic.dv.dims.n <- 4
+
+# 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-------------------------------------------
+
+# start debugging-------------------------------------------
+} else if (input.epic.dv.dims.n == 4) {
+  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-------------------------------------------
+
+} 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
 
 ### Get input data
-input.epic.arr <-
-  ncdf4::ncvar_get(input.epic.fh, varid=input_epic_dv_name, start=input.epic.read.start, count=input.epic.read.count)
+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, I'm caught between IOAPI (which demands tstep=1) and netCDF (which strips tstep=1)
+if ((length(dim(input.epic.dv.arr)) == 3) && (input.epic.dv.dims.n == 4)) {
+  dim(input.epic.dv.arr) <- c(dim(input.epic.dv.arr), 1)
+}
 
 # # 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 emittivities=', this_fn)
 # )
 # #         sum=1.39e+03
 # #   end debugging-----------------------------------------------------
 
+### input TFLAG
+
+## ... 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 datavar=%s',
+    this_fn, input.epic.TFLAG.dims.n, input_epic_TFLAG_name))
+
+} 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 datavar=%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 datavar=%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 datavar=%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 datavar=%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 datavar=%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
+
 # ----------------------------------------------------------------------
 # setup output data container(s)
 # ----------------------------------------------------------------------
 
-# dim(output.epic.arr) <- output.epic.dv.dims
-# gotta create arrays explicitly
-output.epic.arr <- array(NA, output.epic.dv.dims)
-output.epic.write.start <- rep(1, output.epic.dv.dims.n)   # start=(1,1,1,...)
-output.epic.write.count <- output.epic.dv.dims             # count=(n,n,n,...)
+### TFLAG
+
+output.epic.TFLAG.dims <- input.epic.TFLAG.dims # (TSTEP, VAR, DATE-TIME), except ...
+output.epic.TFLAG.dims.n <- length(output.epic.TFLAG.dims)
+output.epic.TFLAG.arr <- input.epic.TFLAG.arr
+output.epic.TFLAG.write.start <- input.epic.TFLAG.read.start
+output.epic.TFLAG.write.count <- input.epic.TFLAG.read.count
+
+# # start debugging-------------------------------------------
+# cat(sprintf('%s: output.epic.TFLAG.dims==', this_fn))
+# print(output.epic.TFLAG.dims)
+# cat(sprintf('%s: output.epic.TFLAG.write.start==', this_fn))
+# print(output.epic.TFLAG.write.start)
+# cat(sprintf('%s: output.epic.TFLAG.write.count==', this_fn))
+# print(output.epic.TFLAG.write.count)
+# #   end debugging-------------------------------------------
+
+### the real datavar
+
+output.epic.dv.arr <- array(NA, output.epic.dv.dims)
+output.epic.dv.write.start <- rep(1, output.epic.dv.dims.n)   # start=(1,1,1,...)
+output.epic.dv.write.count <- output.epic.dv.dims             # count=(n,n,n,...)
 
 # # start debugging-----------------------------------------------------
-# base::cat(base::sprintf('\n%s: output.epic.write.start==', this_fn))
-# print(output.epic.write.start)
-# base::cat(base::sprintf('%s: output.epic.write.count==', this_fn))
-# print(output.epic.write.count)
+# base::cat(base::sprintf('\n%s: output.epic.dv.write.start==', this_fn))
+# print(output.epic.dv.write.start)
+# base::cat(base::sprintf('%s: output.epic.dv.write.count==', this_fn))
+# print(output.epic.dv.write.count)
 # #   end debugging-----------------------------------------------------
 
 # ----------------------------------------------------------------------
 
 ### Process each gridcell. TODO: as dual `apply`?
 
-for (i.col in 1:input.epic.dim.x.n) {
-  for (i.row in 1:input.epic.dim.y.n) {
-
-    input.beld.arr[i.col,i.row,] -> input.beld.vec # vector of crop coverages
-    input.epic.arr[i.col,i.row,] -> input.epic.vec # vector of crop emittivities
-
-    emissions.list <- sum.emissions.for.layers(
-      epic.vec=input.epic.vec,
-      beld.vec=input.beld.vec,
-      col=i.col,
-      row=i.row,
-      record.gridcells=record.gridcells,
-      record.have.epic.not.beld=record.have.epic.not.beld,
-      record.have.beld.not.epic=record.have.beld.not.epic
-    )
-    output.epic.arr[i.col,i.row] <- emissions.list$emissions # scalar product of input vectors
-    # sum.emissions.for.layers "rewrites" statistical/reporting inputs
-    record.gridcells <- emissions.list$gridcells.vec
-    record.have.epic.not.beld <- emissions.list$have.epic.not.beld
-    record.have.beld.not.epic <- emissions.list$have.beld.not.epic
-
-  } # end for rows
-} # end for cols
+for (i.tstep in 1:input.epic.dim.tsteps.n) {
+# i.tstep <- 1
+
+  input.epic.dv.arr.tstep <- input.epic.dv.arr[,,,i.tstep]
+#   for (i.lay in 1:input.epic.dim.layers.n) {
+    for (i.col in 1:input.epic.dim.x.n) {
+      for (i.row in 1:input.epic.dim.y.n) {
+        # both input.beld.vec and input.epic.vec extend along LAY
+        input.beld.vec <- input.beld.arr[i.col,i.row,]       # vector of crop coverages
+        input.epic.vec <- input.epic.dv.arr.tstep[i.col,i.row,] # vector of crop emittivities
+
+        emissions.list <- sum.emissions.for.layers(
+          epic.vec=input.epic.vec,
+          beld.vec=input.beld.vec,
+          col=i.col,
+          row=i.row,
+          record.gridcells=record.gridcells,
+          record.have.epic.not.beld=record.have.epic.not.beld,
+          record.have.beld.not.epic=record.have.beld.not.epic
+        )
+#        output.epic.dv.arr[i.col, i.row, i.lay, i.tstep] <-
+# ASSERT: output has only one layer=LAY!
+        output.epic.dv.arr[i.col, i.row, 1, i.tstep] <-
+          emissions.list$emissions # scalar product of input vectors
+        # sum.emissions.for.layers "rewrites" statistical/reporting inputs
+        record.gridcells <- emissions.list$gridcells.vec
+        record.have.epic.not.beld <- emissions.list$have.epic.not.beld
+        record.have.beld.not.epic <- emissions.list$have.beld.not.epic
+
+      } # end for rows
+    } # end for cols
+#   } # end for layers
+} # end for tsteps
 
 # start debugging-----------------------------------------------------
-base::cat(base::sprintf('\n%s: dim(output.epic.arr)==', this_fn))
-print(dim(output.epic.arr))
+base::cat(base::sprintf('\n%s: dim(output.epic.dv.arr)==', this_fn))
+print(dim(output.epic.dv.arr))
 stats.to.stdout(
-  data=output.epic.arr,
+  data=output.epic.dv.arr,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: EPIC output emissions=', this_fn)
 )
 
 # Very big gap between q3=0.00111 and max=0.283: gotta suspect an outlier/error.
 # What are our results minus that max?
-output.epic.arr.max <- max(output.epic.arr, na.rm=TRUE)
-output.epic.arr.max.coords <- # returns 2D: useful to humans, not for R
-  as.vector(which(output.epic.arr >= output.epic.arr.max, arr.ind = TRUE))
-output.epic.arr.max.index <- # returns 1D: R can use to lookup/change the value
-  which(output.epic.arr >= output.epic.arr.max)
+output.epic.dv.arr.max <- max(output.epic.dv.arr, na.rm=TRUE)
+output.epic.dv.arr.max.coords <- # returns 2D: useful to humans, not for R
+  as.vector(which(output.epic.dv.arr >= output.epic.dv.arr.max, arr.ind = TRUE))
+output.epic.dv.arr.max.index <- # returns 1D: R can use to lookup/change the value
+  which(output.epic.dv.arr >= output.epic.dv.arr.max)
 
 # # start debugging-----------------------------------------------------
-# base::cat(base::sprintf('\n%s: max(output.epic.arr)==%e\n',
-#   this_fn, output.epic.arr.max))
-# base::cat(base::sprintf('\n%s: coords max(output.epic.arr)==', this_fn))
-# print(output.epic.arr.max.coords)
+# base::cat(base::sprintf('\n%s: max(output.epic.dv.arr)==%e\n',
+#   this_fn, output.epic.dv.arr.max))
+# base::cat(base::sprintf('\n%s: coords max(output.epic.dv.arr)==', this_fn))
+# print(output.epic.dv.arr.max.coords)
 
-# # compute_EPIC_emissions.r: max(output.epic.arr)==2.829876e-01
-# # compute_EPIC_emissions.r: coords max(output.epic.arr)==[1]  93 235
+# # compute_EPIC_emissions.r: max(output.epic.dv.arr)==2.829876e-01
+# # compute_EPIC_emissions.r: coords max(output.epic.dv.arr)==[1]  93 235
 # #   end debugging-----------------------------------------------------
 
 # results minus that max
-output.epic.arr.minus.max <- output.epic.arr
-output.epic.arr.minus.max[output.epic.arr.max.index] <- NA
+output.epic.dv.arr.minus.max <- output.epic.dv.arr
+output.epic.dv.arr.minus.max[output.epic.dv.arr.max.index] <- NA
 
 # start debugging-----------------------------------------------------
 stats.to.stdout(
-  data=output.epic.arr.minus.max,
+  data=output.epic.dv.arr.minus.max,
   sig.digs=sigdigs,
   title=base::sprintf('\n%s: EPIC output emissions minus (spurious?) max=%.3g:',
-    this_fn, output.epic.arr.max)
+    this_fn, output.epic.dv.arr.max)
 )
 
 # compute_EPIC_emissions.r: EPIC output emissions minus (spurious?) max=0.283:
 # start debugging-----------------------------------------------------
 base::cat(base::sprintf(
   '\n%s: WARNING: removing presumed-spurious max=%e from output\n',
-  this_fn, output.epic.arr.max))
+  this_fn, output.epic.dv.arr.max))
 #   end debugging-----------------------------------------------------
 
-output.epic.arr <- output.epic.arr.minus.max
+output.epic.dv.arr <- output.epic.dv.arr.minus.max
+
+# ----------------------------------------------------------------------
+# 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.epic.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
+cat(sprintf('\n%s: 1. define output dimensions\n', this_fn)) # debugging
 
-## instead, copy input dimensions? No: above have type=list, not type=ncdim4 :-(
+## instead, copy input dimensions? I can't see how to :-(
 
-# input.epic.dim.x.vals <- input.epic.dim.x$vals
-# input.epic.dim.y.vals <- input.epic.dim.y$vals
-## To overlay plots with "the usual" map, need dimension centered and units=m
+## 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 <- 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):
+# > [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 <- 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.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.vals <- input.epic.dim.vars$vals
+# input.epic.dim.vars.units <- input.epic.dim.vars$units
+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.epic.dim.x <- ncdim_def(
   name=output_epic_dim_x_name,
-#   units='',
   units=output_epic_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_epic_dim_x_long_name)
 
 output.epic.dim.y <- ncdim_def(
   name=output_epic_dim_y_name,
-#   units='',
   units=output_epic_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_epic_dim_y_long_name)
 
+# output=input
+output.epic.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.epic.dim.layers <- ncdim_def(
+  name=input_epic_dim_layers_name,
+  units=input.epic.dim.layers.units,
+#  vals=input.epic.dim.layers.vals, # 1:42, not what output wants
+#  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.epic.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.epic.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
+
 ### 2. define output datavar(s)
-# cat(sprintf('\n%s: 2. define output datavar(s)\n', this_fn)) # debugging
+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.epic.TFLAG <- ncdf4::ncvar_def(
+  name=output_epic_TFLAG_name,
+  units=input.epic.TFLAG.units,
+# (TSTEP, VAR, DATE-TIME) is how it looks in `ncdump`, but ncdf4 wants ...
+  dim=list(output.epic.dim.datetimes, output.epic.dim.vars, output.epic.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.epic.dv <- ncdf4::ncvar_def(
   name=output_epic_dv_name,
   units=output_epic_dv_units,
-#  dim=list(output.epic.dim.y, output.epic.dim.x),
-# (ROW, COL) is how it looks in `ncdump`, but internally we want ...
-  dim=list(output.epic.dim.x, output.epic.dim.y),
+# (TSTEP, LAY, ROW, COL) is how it looks in `ncdump`, but ncdf4 wants ...
+  dim=list(output.epic.dim.x, output.epic.dim.y, output.epic.dim.layers, output.epic.dim.tsteps),
   missval=input.epic.dv.mv,
   longname=output_epic_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
+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.epic.fh <- ncdf4::nc_create(
   filename=output_epic_fp,
-  vars=output.epic.dv,               # can also be vector or list
+#   vars=output.epic.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.epic.TFLAG, output.epic.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
+cat(sprintf('\n%s: 3.1. create output file/global attributes\n', this_fn)) # debugging
 
 ## copy attributes from EPIC input, then overwrite {NLAYS, VAR-LIST, FILEDESC}
 
 ## After much thrashing, and some guidance from http://www.r-bloggers.com/s4-classes-in-r-printing-function-definition-and-getting-help/
 ## I discover howto get the function definition: `package:::methodname`, e.g., `ncdf4:::print.ncdf4`
 ## which shows the secret ncdf4 syntax for getting/walking all the file/global attributes:
-## > atts <- ncatt_get(nc, 0)
+## > atts <- ncatt_get(nc, 0) # where `nc`==filehandle
 ## > natts <- length(atts)
 ## > if (natts > 0) {
 ## >   print(paste("    ", natts, " global attributes:", sep = ""))
 } # 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
+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
+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.epic.fh == %s\n', output_epic_fp))
 # base::cat(base::sprintf('  varid=%s\n', output_epic_dv_name))
-# base::cat(base::sprintf('  vals=output.epic.arr == '))
-# print(dim(output.epic.arr))
+# base::cat(base::sprintf('  vals=output.epic.dv.arr == '))
+# print(dim(output.epic.dv.arr))
 # #   end debugging-----------------------------------------------------
 
 ncdf4::ncvar_put(
   nc=output.epic.fh,
+#   varid=output_epic_TFLAG_name,
+  varid=output.epic.TFLAG,
+  vals=output.epic.TFLAG.arr,
+  start=output.epic.TFLAG.write.start,
+  count=output.epic.TFLAG.write.count,
+# # to debug, toggle verbose
+#   verbose=TRUE)
+  verbose=FALSE)
+
+ncdf4::ncvar_put(
+  nc=output.epic.fh,
 #  varid=output.epic.dv,
   varid=output_epic_dv_name,
-  vals=output.epic.arr,
-  start=output.epic.write.start,
-  count=output.epic.write.count,
+  vals=output.epic.dv.arr,
+  start=output.epic.dv.write.start,
+  count=output.epic.dv.write.count,
 # # to debug, toggle verbose
 #   verbose=TRUE)
   verbose=FALSE)
 
 ## 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
+cat(sprintf('\n%s: 4.2. write datavar metadata (i.e., additional attributes)\n', this_fn)) # debugging
+
+ncdf4::ncatt_put(
+  nc=output.epic.fh,
+  varid=output.epic.TFLAG,
+  attname='var_desc',
+  attval=output.epic.TFLAG.var_desc,
+  prec='text',
+  definemode=FALSE,
+# # to debug, toggle verbose
+#   verbose=TRUE)
+  verbose=FALSE)
 
 ncdf4::ncatt_put(
   nc=output.epic.fh,
   layers.n=input.epic.dim.layers.n         # total number of layers
 )
 
+#----------------------------------------------------------------------
+# verify VERDI compatibility? output=EPIC emissions
+#----------------------------------------------------------------------
+
+  ### TODO: integrate VERDI into bash_utilities.sh
+  ### plot our datavar of interest with VERDI
+  output_epic_dv_name_VERDI <- sprintf('%s[1]', output_epic_dv_name)
+  ## VERDI really wants an FQP
+  output_epic_fp_VERDI <- system(sprintf('readlink -f %s', output_epic_fp), intern=TRUE)
+  output_epic_dv_VERDI_cmd <-
+    sprintf('verdi -f %s -s %s -gtype tile &',
+      output_epic_fp_VERDI, output_epic_dv_name_VERDI)
+  cat(sprintf("\n%s: debugging IOAPI: about to run '%s'\n", this_fn, output_epic_dv_VERDI_cmd))
+#   cat(sprintf("\n%s: not debugging IOAPI: skipping '%s'\n", this_fn, output_epic_dv_VERDI_cmd))
+  system(output_epic_dv_VERDI_cmd)
+
 # ----------------------------------------------------------------------
 # visualize output emissions
 # ----------------------------------------------------------------------

distill_EPIC_emittivities.ncl

-;!/usr/bin/env ncl ; requires version >= ???
+;!/usr/bin/env ncl ; requires version >= 6.0.0 for default_fillvalue(...)
 ;----------------------------------------------------------------------
 ;;; TODO: describe me!
 
 
 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 getFillValue, isMonotonic
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; for getFillValue, isMonotonic, copy_VarMeta
 
-; load "$WORK_DIR/get_daily_emis_fp.ncl" ; for that function
-load "$WORK_DIR/string.ncl" ; for pad_blanks_to_length
-; load "$WORK_DIR/summarize.ncl" ; for
-load "$WORK_DIR/time.ncl" ; for seconds_in_year
+load "$STRING_FUNCS_FP" ; for pad_blanks_to_length
+load "$TIME_FUNCS_FP" ; for seconds_in_year
 
 ;----------------------------------------------------------------------
 ; constants
   ret_vec = new(dimsizes(raw_vec), raw_type, getFillValue(raw_array))
   ret_vec(dense_vec_indices) = demono_vec
   ; ... and transform that to shape of input
-
   ret_array = onedtond(ret_vec, dimsizes(raw_array))
 
 ;   ; start debugging---------------------------------------------------
 ; load only species=N2O
 ;----------------------------------------------------------------------
 
-  ;;; raw* (other than file and filehandle) have only one species (N2O), possibly monotonicity, and flux-rate units
+  ;;; raw* may have multiple species (other than N2O), possibly monotonicity, and definitely flux-rate units
   raw_emitt_fh = addfile(raw_emitt_fp, "r")
   raw_emitt_dv = raw_emitt_fh->$raw_emitt_dv_name$
   raw_emitt_dv_dims_sizes = dimsizes(raw_emitt_dv)
   raw_emitt_dv_dims_n = dimsizes(raw_emitt_dv_dims_sizes) ; == |dims|
 ;   raw_emitt_dv_dims_names = (/ "TSTEP", "LAY", "ROW", "COL" /) ; TODO: get programmatically
   raw_emitt_dv_dims_names = getvardims(raw_emitt_dv)
+  raw_emitt_dv_tsteps_n = raw_emitt_dv_dims_sizes(0)
   raw_emitt_dv_layers_n = raw_emitt_dv_dims_sizes(1) ; LAY == crop
   raw_emitt_dv_rows_n = raw_emitt_dv_dims_sizes(2)
   raw_emitt_dv_cols_n = raw_emitt_dv_dims_sizes(3)
   raw_emitt_dv_type = typeof(raw_emitt_dv)     ; type as string
 ;   raw_emitt_dv_units = raw_emitt_dv@units
 
-  ;;; distilled* have only N2O, definitely no monotonicity, and flux-rate units
+;----------------------------------------------------------------------
+; copy/mod IOAPI/fake datavar=TFLAG
+;----------------------------------------------------------------------
+
+  ;;; For IOAPI compatibility (mostly so EJC can load this into VERDI), we need IOAPI's fake "datavar" name=TFLAG
+  raw_emitt_TFLAG = raw_emitt_fh->TFLAG
+  raw_emitt_TFLAG_type = typeof(raw_emitt_TFLAG)     ; type as string
+  raw_emitt_TFLAG_mv = getFillValue(raw_emitt_TFLAG) ; missing value
+
+;;; Major aside: IOAPI does not apparently typically set attr={_FillValue, missing_value} on datavar=TFLAG, which makes raw_emitt_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 raw_emitt_TFLAG_mv, and
+;;; TODO: raise issue on netcdfgroup@unidata.ucar.edu: NCL OK with putting datavar with missing _FillValue, ncdf4 is not
+
+  if (ismissing(raw_emitt_TFLAG_mv)) then
+    raw_emitt_TFLAG_mv = default_fillvalue(raw_emitt_TFLAG_type)
+  end if
+  raw_emitt_TFLAG_dims_sizes = dimsizes(raw_emitt_TFLAG)
+  raw_emitt_TFLAG_dims_n = dimsizes(raw_emitt_TFLAG_dims_sizes) ; == |dims|
+  raw_emitt_TFLAG_dims_names = getvardims(raw_emitt_TFLAG) ; TSTEP, VAR, DATE-TIME
+
+;   ; start debugging----------------------------------------------------
+;   print("printVarSummary(raw_emitt_TFLAG)==")
+;   printVarSummary(raw_emitt_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+  ;;; output should have
+  ;;; * same |TSTEP|
+  pure_emitt_TFLAG_tsteps_n = raw_emitt_TFLAG_dims_sizes(0)
+  ;;; * |VAR| = 1 (only one "real" datavar=N2O)
+  pure_emitt_TFLAG_vars_n = 1
+  ;;; * same |DATE-TIME|
+  pure_emitt_TFLAG_datetimes_n = raw_emitt_TFLAG_dims_sizes(2)
+  pure_emitt_TFLAG_dims_sizes = \
+    (/ pure_emitt_TFLAG_tsteps_n, pure_emitt_TFLAG_vars_n, pure_emitt_TFLAG_datetimes_n /)
+
+  ;;; make the "datavar", and ...
+  pure_emitt_TFLAG = \
+    new(pure_emitt_TFLAG_dims_sizes, raw_emitt_TFLAG_type, raw_emitt_TFLAG_mv)
+  ;;; ... copy its dimensions and attributes
+  copy_VarMeta(raw_emitt_TFLAG, pure_emitt_TFLAG)
+
+;   ; start debugging----------------------------------------------------
+;   print(this_fn+": pre-population: printVarSummary(pure_emitt_TFLAG)==")
+;   printVarSummary(pure_emitt_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+  ;;; iterate (TSTEP, VAR, DATE-TIME)
+
+;   do i_tstep = 0, 0
+  do i_tstep = 0, pure_emitt_TFLAG_tsteps_n-1
+;     do i_var = 0, 0 ; is all we're gonna do, but ...
+    do i_var = 0, pure_emitt_TFLAG_vars_n-1
+      ;; Note that TFLAG values are in fact VAR-independent: see https://github.com/TomRoche/ioapi "Step 4"
+      ;; so we only need to copy one
+      pure_emitt_TFLAG(i_tstep, i_var, :) = raw_emitt_TFLAG(i_tstep, i_var, :)
+    end do ; i_var = 0, raw_emitt_TFLAG_vars_n-1
+  end do ; i_tstep = 0, raw_emitt_TFLAG_tsteps_n-1
+
+;   ; start debugging----------------------------------------------------
+;   print(this_fn+": post-population: print(pure_emitt_TFLAG)==")
+;   print(pure_emitt_TFLAG)
+;   ;   end debugging----------------------------------------------------
+
+;----------------------------------------------------------------------
+; load species=N2O && TFLAG for IOAPI/VERDI
+;----------------------------------------------------------------------
+
+  ;;; distilled* have only one species=N2O (and the IOAPI fake datavar=TFLAG), definitely no monotonicity, and flux-rate units
   ;;; raw_emitt_dv has dimensions=(TSTEP, LAY, ROW, COL)
   ;;; * but there's only one TSTEP (TODO: test that!)
-  ;;; so ignore TSTEP
-  dist_emitt_dv_dims_sizes = raw_emitt_dv_dims_sizes(1:)
-  dist_emitt_dv_dims_names = raw_emitt_dv_dims_names(1:)
-  dist_emitt_dv_dims_n = raw_emitt_dv_dims_n - 1
-  dist_emitt_dv_layers_n = dist_emitt_dv_dims_sizes(0)
-  dist_emitt_dv_rows_n = dist_emitt_dv_dims_sizes(1)
-  dist_emitt_dv_cols_n = dist_emitt_dv_dims_sizes(2)
+  ;;; so ignore TSTEP? not for IOAPI/VERDI compatibility
+;   dist_emitt_dv_dims_sizes = raw_emitt_dv_dims_sizes(1:)
+;   dist_emitt_dv_dims_names = raw_emitt_dv_dims_names(1:)
+;   dist_emitt_dv_dims_n = raw_emitt_dv_dims_n - 1
+;   dist_emitt_dv_layers_n = dist_emitt_dv_dims_sizes(0)
+;   dist_emitt_dv_rows_n = dist_emitt_dv_dims_sizes(1)
+;   dist_emitt_dv_cols_n = dist_emitt_dv_dims_sizes(2)
+
+; striving for IOAPI/VERDI compliance
+  dist_emitt_dv_dims_sizes = raw_emitt_dv_dims_sizes
+  dist_emitt_dv_dims_names = raw_emitt_dv_dims_names
+  dist_emitt_dv_dims_n = raw_emitt_dv_dims_n
+  dist_emitt_dv_tsteps_n = raw_emitt_dv_tsteps_n
+  dist_emitt_dv_layers_n = raw_emitt_dv_layers_n ; yes! not 1: that happens when we compute emissions, not now.
+  dist_emitt_dv_rows_n = raw_emitt_dv_rows_n
+  dist_emitt_dv_cols_n = raw_emitt_dv_cols_n
   dist_emitt_arr = \
     new(dist_emitt_dv_dims_sizes, raw_emitt_dv_type, raw_emitt_dv_mv)
 
 ;   printVarSummary(areas_arr)
 ;   ;   end debugging---------------------------------------------------
 
-  ;;; pure* have only one species (N2O), no monotonicity, and mass-rate units
+  ;;; pure* have only one species (N2O) (and the IOAPI fake datavar=TFLAG), no monotonicity, and mass-rate units
   pure_emitt_dv_dims_sizes = dist_emitt_dv_dims_sizes
   pure_emitt_dv_dims_names = dist_emitt_dv_dims_names
   pure_emitt_dv_dims_n = dist_emitt_dv_dims_n
   pure_emitt_arr = \
     new(pure_emitt_dv_dims_sizes, raw_emitt_dv_type, raw_emitt_dv_mv)
 
+  ;;; get dimensions for the whole file, since != dimensions for either datavar
+  pure_emitt_fh_dims_names = getvardims(raw_emitt_fh)
+  pure_emitt_fh_dims_n = dimsizes(pure_emitt_fh_dims_names)
+  ;;; more complex
+  pure_emitt_fh_dims_sizes = new(pure_emitt_fh_dims_n, integer, "No_FillValue")
+  do i_dim = 0, pure_emitt_fh_dims_n - 1
+    i_name = pure_emitt_fh_dims_names(i_dim)
+    if      (i_name .eq. "COL") then
+      pure_emitt_fh_dims_sizes(i_dim) = dist_emitt_dv_cols_n
+    else if (i_name .eq. "DATE-TIME") then
+      pure_emitt_fh_dims_sizes(i_dim) = pure_emitt_TFLAG_datetimes_n
+    else if (i_name .eq. "LAY") then
+      pure_emitt_fh_dims_sizes(i_dim) = dist_emitt_dv_layers_n
+    else if (i_name .eq. "ROW") then
+      pure_emitt_fh_dims_sizes(i_dim) = dist_emitt_dv_rows_n
+    else if (i_name .eq. "VAR") then
+      pure_emitt_fh_dims_sizes(i_dim) = pure_emitt_TFLAG_vars_n
+    else if (i_name .eq. "TSTEP") then
+      pure_emitt_fh_dims_sizes(i_dim) = dist_emitt_dv_tsteps_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 ; i_dim = 0, pure_emitt_fh_dims_n - 1
+
+  ; start debugging---------------------------------------------------
+  print("print(pure_emitt_fh_dims_names)==")
+  print(pure_emitt_fh_dims_names)
+  print("print(pure_emitt_fh_dims_sizes)==")
+  print(pure_emitt_fh_dims_sizes)
+  ;   end debugging---------------------------------------------------
+
+;----------------------------------------------------------------------
+; verify VERDI compatibility? raw input
+;----------------------------------------------------------------------
+
+  ;;; TODO: integrate VERDI into bash_utilities.sh
+  ;;; plot our datavar of interest with VERDI
+  raw_emitt_dv_name_VERDI = raw_emitt_dv_name+"[1]"
+  ;; VERDI really wants an FQP
+  raw_emitt_fp_VERDI = systemfunc("readlink -f "+raw_emitt_fp)
+  raw_emitt_dv_VERDI_cmd = "verdi -f "+raw_emitt_fp_VERDI+" -s "+raw_emitt_dv_name_VERDI+" -gtype tile &"
+;   print(str_get_nl() + \
+;     this_fn+": debugging IOAPI: about to run '"+raw_emitt_dv_VERDI_cmd+"'" + str_get_nl())
+;   system(raw_emitt_dv_VERDI_cmd)
+  print(str_get_nl() + \
+    this_fn+": not debugging IOAPI: skipping '"+raw_emitt_dv_VERDI_cmd+"'" + str_get_nl())
+
 ;----------------------------------------------------------------------
 ; demonotonicize emittivities if necessary
 ;----------------------------------------------------------------------
 
-;;; iterate (layers, rows, cols)
-;   do i_layer = 0, 0
-  do i_layer = 0, raw_emitt_dv_layers_n-1
-    layer_raw = raw_emitt_dv(0,i_layer,:,:)
-    layer_raw_monotonicity = monotonicity(layer_raw)
-    if (layer_raw_monotonicity .eq. 0) then
-      ; not monotonic? take as-is
-      dist_emitt_arr(i_layer,:,:) = layer_raw
-    else if (layer_raw_monotonicity .eq. 1) then
-      ; increasing monotonicity? remove it
-;       print(this_fn+": DEBUGGING: layer="+i_layer+" monotonic increasing")
-      dist_emitt_arr(i_layer,:,:) = (/ remove_increasing_monotonicity(layer_raw) /)
-    else
-      ; TODO: throw error
-      print(str_get_nl() + \
-        this_fn+": ERROR: panic: decreasing monotonicity" + \
-        str_get_nl())
-      return ; TODO: trainwreck
-    end if
-    end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
-  end do ; i_layer = 0, raw_emitt_dv_layers_n-1
+;;; iterate (tsteps, layers, rows, cols)
+
+;   do i_tstep = 0, 0
+  do i_tstep = 0, raw_emitt_dv_tsteps_n-1
+;     do i_layer = 0, 0
+    do i_layer = 0, raw_emitt_dv_layers_n-1
+      layer_raw = raw_emitt_dv(i_tstep, i_layer, :,:) ; all rows, cols
+      layer_raw_monotonicity = monotonicity(layer_raw)
+      if (layer_raw_monotonicity .eq. 0) then
+        ; not monotonic? take as-is
+        dist_emitt_arr(i_tstep,i_layer,:,:) = layer_raw
+      else if (layer_raw_monotonicity .eq. 1) then
+        ; increasing monotonicity? remove it
+;         print(this_fn+": DEBUGGING: layer="+i_layer+" monotonic increasing")
+        dist_emitt_arr(i_tstep,i_layer,:,:) = (/ remove_increasing_monotonicity(layer_raw) /)
+      else
+        ; TODO: throw error
+        print(str_get_nl() + \
+          this_fn+": ERROR: panic: decreasing monotonicity" + \
+          str_get_nl())
+        return ; TODO: trainwreck
+      end if
+      end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+    end do ; i_layer = 0, raw_emitt_dv_layers_n-1
+  end do ; i_tstep = 0, raw_emitt_dv_tsteps_n-1  
 
 ;   ; start debugging----------------------------------------------------
 ;   print("printVarSummary(dist_emitt_arr)==")
   ;;; ------ = ----- * --- * -------- * -- * ----- * ------ 
   ;;;      s   ha yr   m^2   gridcell   s    kgN2O   gN2O    
 
-;;; iterate (layers, rows, cols)
-;   do i_layer = 0, 0
-  do i_layer = 0, dist_emitt_dv_layers_n-1
-    layer_dist = dist_emitt_arr(i_layer,:,:)
-    pure_emitt_arr(i_layer,:,:) = \
-      layer_dist         * \ ; kgN20/ha/yr
-      1e-4               * \ ; ha/m^2
-      areas_arr          * \ ; m^2/gridcell
-      years_per_second   * \ ; yr/s
-      1e3                * \ ; gN2O/kgN2O
-      molN2O_per_gN2O        ; molN2O/gN2O
-  end do ; i_layer = 0, dist_emitt_dv_layers_n-1
+;;; iterate (tsteps, layers, rows, cols)
+
+;   do i_tstep = 0, 0
+  do i_tstep = 0, raw_emitt_dv_tsteps_n-1
+;     do i_layer = 0, 0
+    do i_layer = 0, dist_emitt_dv_layers_n-1
+      layer_dist = dist_emitt_arr(i_tstep,i_layer,:,:)
+      pure_emitt_arr(i_tstep,i_layer,:,:) = \
+        layer_dist         * \ ; kgN20/ha/yr
+        1e-4               * \ ; ha/m^2
+        areas_arr          * \ ; m^2/gridcell
+        years_per_second   * \ ; yr/s
+        1e3                * \ ; gN2O/kgN2O
+        molN2O_per_gN2O        ; molN2O/gN2O
+    end do ; i_layer = 0, dist_emitt_dv_layers_n-1
+  end do ; i_tstep = 0, raw_emitt_dv_tsteps_n-1  
 
 ;   ; start debugging----------------------------------------------------
 ;   print("printVarSummary(pure_emitt_arr)==")
   if (.not. all(ismissing(raw_emitt_fa_names))) then
     do i=0, dimsizes(raw_emitt_fa_names)-1
       fa_name = raw_emitt_fa_names(i)
-      if      (fa_name .eq. "NVARS") then
+;       if      (fa_name .eq. "NLAYS") then
+;         pure_emitt_fh@$fa_name$ = 1 ; |data variables|
+;       else if (fa_name .eq. "NVARS") then
+; NO! NLAYS -> 1 when we do emissions, not now
+      if (fa_name .eq. "NVARS") then
         pure_emitt_fh@$fa_name$ = 1 ; |data variables|
       else if (fa_name .eq. "VAR-LIST") then
         ; IOAPI-style
         pure_emitt_fh@$fa_name$ = pad_blanks_to_length(pure_emitt_dv_name, 16)
       else
         pure_emitt_fh@$fa_name$ = raw_emitt_fh@$fa_name$
+;       end if
       end if
       end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
     end do
   pure_emitt_fh@HISTORY = pure_emitt_attr_history
 
   ;;; 3. Define {dimensions, coordinate variables} and their dimensionality
-  ;;; don't copy all from input: we only want those needed for pure_emitt_arr
+
+  ;;; don't copy all from input? we only want those needed for pure_emitt_arr?
 
   ;;; 3.1. Define output dimensions.
   ;; Should any dimensions be unlimited? Here, no.
-;  pure_emitt_dv_dims_unlim = (/ False * pure_emitt_dv_dims_n /) ; too simple :-(
-  pure_emitt_dv_dims_unlim = new(pure_emitt_dv_dims_n, logical, "No_FillValue")
-  do i=0 , pure_emitt_dv_dims_n-1
-    pure_emitt_dv_dims_unlim(i) = False
+  pure_emitt_fh_dims_unlim = new(pure_emitt_fh_dims_n, logical, "No_FillValue")
+  do i=0 , pure_emitt_fh_dims_n-1
+    pure_emitt_fh_dims_unlim(i) = False
   end do
-  filedimdef(pure_emitt_fh, pure_emitt_dv_dims_names, pure_emitt_dv_dims_sizes, pure_emitt_dv_dims_unlim)
-
-  ;;; 3.2. Create output coordinate variables (if they exist in the input).
-  ;;; procedure adapted from contributed.ncl::copy_VarCoords
-  if (.not. all(ismissing(raw_emitt_dv_dims_names))) then
-; Skip this for our case: we explicitly want |output dimensions| < |input dimensions|
-;     if (all(raw_emitt_dv_dims_sizes(:pure_emitt_dv_dims_n-1) .eq. pure_emitt_dv_dims_sizes)) then
-    do i=0, pure_emitt_dv_dims_n-1
-; iterate over the output dimension names, which omit the first input dimension (name)
-      if (.not. ismissing(pure_emitt_dv_dims_names(i))) then
-        pure_emitt_arr!i = raw_emitt_dv!(i+1)
-        if (iscoord(raw_emitt_dv, raw_emitt_dv!(i+1)))
-          pure_emitt_arr&$pure_emitt_arr!i$ = raw_emitt_dv&$raw_emitt_dv!i$
-        end if
-      end if
-    end do
-;   else
-;     print(this_fn+": ERROR: dimension sizes do not match: input="+raw_emitt_dv_dims_sizes+", output="+pure_emitt_dv_dims_sizes)
-;   end if ; checking in and out dim sizes
-  end if ; .not. all(ismissing(raw_emitt_dv_dims_names))
+;  filedimdef(pure_emitt_fh, pure_emitt_dv_dims_names, pure_emitt_dv_dims_sizes, pure_emitt_dv_dims_unlim)
+; NO! we want all the dims for the *file*
+;   filedimdef(pure_emitt_fh, pure_emitt_fh_dims_names, pure_emitt_dv_dims_sizes, pure_emitt_fh_dims_unlim)
+; this fails with
+; > ncdimdef: ncid 196608: NetCDF: Invalid dimension size
+; ???
+  filedimdef(pure_emitt_fh,         \
+    (/ pure_emitt_fh_dims_names /), \
+    (/ pure_emitt_fh_dims_sizes /), \
+    (/ pure_emitt_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(raw_emitt_dv_dims_names))) then
+; ; Skip this for our case: we explicitly want |output dimensions| < |input dimensions|
+; ;     if (all(raw_emitt_dv_dims_sizes(:pure_emitt_dv_dims_n-1) .eq. pure_emitt_dv_dims_sizes)) then
+;     do i=0 /), pure_emitt_dv_dims_n-1
+; ; iterate over the output dimension names, which omit the first input dimension (name)
+;       if (.not. ismissing(pure_emitt_dv_dims_names(i))) then
+;         pure_emitt_arr!i = raw_emitt_dv!(i+1)
+;         if (iscoord(raw_emitt_dv, raw_emitt_dv!(i+1)))
+;           pure_emitt_arr&$pure_emitt_arr!i$ = raw_emitt_dv&$raw_emitt_dv!i$
+;         end if
+;       end if
+;     end do
+; ;   else
+; ;     print(this_fn+": ERROR: dimension sizes do not match: input="+raw_emitt_dv_dims_sizes+", output="+pure_emitt_dv_dims_sizes)
+; ;   end if ; checking in and out dim sizes
+;   end if ; .not. all(ismissing(raw_emitt_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(pure_emitt_fh, "TFLAG", raw_emitt_TFLAG_type, raw_emitt_TFLAG_dims_names)
+  ; now do the real data
   filevardef(pure_emitt_fh, pure_emitt_dv_name, raw_emitt_dv_type, pure_emitt_dv_dims_names)
 
   ;; 4.2. Define output variable(s) attributes
 
   ; 4.2.1. Copy output variable(s) attributes
+  filevarattdef(pure_emitt_fh, "TFLAG", pure_emitt_TFLAG)
   filevarattdef(pure_emitt_fh, pure_emitt_dv_name, pure_emitt_arr)
 
-  ;; 4.2.2. Overwrite our own output variable(s) attributes
+  ;; 4.2.2. Overwrite our "real" output datavar(s) attributes
   pure_emitt_fh->$pure_emitt_dv_name$@units = pure_emitt_dv_units
   pure_emitt_fh->$pure_emitt_dv_name$@long_name = pure_emitt_dv_long_name
   pure_emitt_fh->$pure_emitt_dv_name$@var_desc = pure_emitt_dv_var_desc
   setfileoption(pure_emitt_fh, "DefineMode", False)
 
   ;;; 5. Write data to output variable(s)
+  pure_emitt_fh->TFLAG = (/ pure_emitt_TFLAG /)
   pure_emitt_fh->$pure_emitt_dv_name$ = (/ pure_emitt_arr /)
 
   ;;; 6. Close output filehandle to flush file: see http://www.ncl.ucar.edu/Support/talk_archives/2012/2196.html
   delete(pure_emitt_fh)
 
+;----------------------------------------------------------------------
+; verify VERDI compatibility? "pure" output
+;----------------------------------------------------------------------
+
+  ;;; TODO: integrate VERDI into bash_utilities.sh
+  ;;; plot our datavar of interest with VERDI
+  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"
+
+;   print(str_get_nl() + \
+;     this_fn+"; debugging IOAPI: about to run '"+pure_emitt_dv_VERDI_cmd+"'" + str_get_nl())
+;   system(pure_emitt_dv_VERDI_cmd)
+  print(str_get_nl() + \
+    this_fn+"; not debugging IOAPI: skipping '"+pure_emitt_dv_VERDI_cmd+"'" + str_get_nl())
+
 end ; distill_EPIC_emittivities.ncl
 out_datavar_long_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DV_LONG_NAME')
 out_datavar_dim_x_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_X_NAME')
 out_datavar_dim_y_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_Y_NAME')
-out_datavar_dim_z_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_TIME_NAME')
-out_datavar_dim_z_units <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_TIME_UNITS')
+out_datavar_dim_z_name <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_NAME')
+out_datavar_dim_z_units <- Sys.getenv('AQMEII_EDGAR_N2O_REGRID_DIM_TSTEP_UNITS')
 
 ### plotting regridded EDGAR data