Commits

Tom Roche committed 42d3e56

fixed row reversal, more thorough output visualization, improved IOAPI

retemp.ncl
* makes fully IOAPI-zed output, which it (currently commented out) views with VERDI
to ensure proper grid orientation: see row-reversal GEIA problem, corrected similarly.
* skips on finding monthly.ncf output

Also:
* developed new regrid_utils functionality, committed there
* removed redundant debugging from check_conservation.ncl

Tested on terrae/EMVL from uber_driver.sh::file in fresh terminal, will repeat with uber_driver.sh::repo soonest

TODO:
* (à la GEIA) remove unnecessary/problematic use of TEMPLATE_EXTENT*, use only TEMPLATE_IOAPI*
* move all these TODOs to https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na/issues
* port *_driver improvements to
** other *_driver: EDGAR_driver, GFED_driver
** other uber_driver: AQMEII::uber_driver, GFED::uber_driver
* 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.*(...)
* fully document platform versions (e.g., linux, compilers, bash, 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 (3)

 # where this is hosted
 PROJECT_WEBSITE='https://bitbucket.org/tlroche/clm_cn_global_to_aqmeii-na'
 
+### AQMEII-NA constants
+
+## datavar 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})"
+# IOAPI pads varattr=units to length=16 with trailing spaces
+AQMEIINA_DV_UNITS="$(printf '%-16s' 'moles/s')"
+# IOAPI pads varattr=var_desc to length=80 with trailing spaces
+# 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_TFLAG_UNITS="$(printf '%-16s' '<YYYYDDD,HHMMSS>')"
+# IOAPI pads varattr=var_desc to length=80 with trailing spaces
+AQMEIINA_TFLAG_VAR_DESC="$(printf '%-80s' 'Timestep-valid flags:  (1) YYYYDDD or (2) HHMMSS')"
+
+## dimensions and their attributes
+AQMEIINA_DIM_LAYER_NAME='LAY'
+AQMEIINA_DIM_LAYER_LONG_NAME='index of layers above surface'
+AQMEIINA_DIM_LAYER_UNITS='unitless'
+
+AQMEIINA_DIM_TSTEP_NAME='TSTEP'
+AQMEIINA_DIM_TSTEP_UNITS='' # they vary!
+AQMEIINA_DIM_TSTEP_LONG_NAME='timestep'
+
+# dimensions we don't really need, but VERDI/IOAPI/TFLAG does
+AQMEIINA_DIM_DATETIME_NAME='DATE-TIME'
+AQMEIINA_DIM_VAR_NAME='VAR'
+
+AQMEIINA_DIM_X_N='459'
+AQMEIINA_DIM_X_NAME='COL'
+# AQMEIINA_DIM_X_UNITS='unitless' # my invention
+AQMEIINA_DIM_X_UNITS='m'
+# AQMEIINA_DIM_X_LONG_NAME='Fortran-style index to grid columns from lower-left origin' # my invention (TODO: CHECK it's not from top)
+AQMEIINA_DIM_X_LONG_NAME='grid-center offset from center of projection'
+
+AQMEIINA_DIM_Y_N='299'
+AQMEIINA_DIM_Y_NAME='ROW'
+# AQMEIINA_DIM_Y_UNITS='unitless' # my invention
+AQMEIINA_DIM_Y_UNITS='m'
+# 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}"
+
 ### for visualization (generally)
 export OUTPUT_SIGNIFICANT_DIGITS='3' # see conservation report below
 # PROJ.4 string for unprojected data
 # no, allow user to override repo code: see `get_bash_utils`
 export BASH_UTILS_FP="${WORK_DIR}/${BASH_UTILS_FN}"
 
+export FILEPATH_FUNCS_FN='get_filepath_from_template.ncl'
+export FILEPATH_FUNCS_FP="${WORK_DIR}/${FILEPATH_FUNCS_FN}"
+
 export IOAPI_FUNCS_FN='IOAPI.ncl'
 export IOAPI_FUNCS_FP="${WORK_DIR}/${IOAPI_FUNCS_FN}"
 
 # export TEMPLATE_INPUT_BAND='1' # `ncks` makes dim=TSTEP first
 export TEMPLATE_DATAVAR_NAME='emi_n2o'
 
+### Template for both `raster` extents and IOAPI-writing
+
+# This file should have all required IOAPI metadata (including the "fake datavar"=TFLAG) for a single datavar=N2O.
+# We'll use it to output IOAPI, so we can check grid alignment in VERDI.
+# It must also have the correct AQMEII-NA extents, to use in `raster` regridding.
+TEMPLATE_IOAPI_GZ_URI='https://bitbucket.org/tlroche/aqmeii_ag_soil/downloads/emis_mole_N2O_2008_12US1_cmaq_cb05_soa_2008ab_08c.ncf.gz' # final yearly output from AQMEII_ag_soil, VERDI-viewable
+# formerly also used
+# TEMPLATE_EXTENT_GZ_URI='https://bitbucket.org/tlroche/geia_regrid/downloads/emis_mole_all_20080101_12US1_cmaq_cb05_soa_2008ab_08c.EXTENTS_INPUT.nc.gz'
+# et al
+TEMPLATE_IOAPI_GZ_FN="$(basename ${TEMPLATE_IOAPI_GZ_URI})"
+TEMPLATE_IOAPI_GZ_FP="${WORK_DIR}/${TEMPLATE_IOAPI_GZ_FN}"
+# TEMPLATE_IOAPI_ROOT="${TEMPLATE_IOAPI_GZ_FN%.*}" # everything left of the LAST dot
+# TEMPLATE_IOAPI_FN="${TEMPLATE_IOAPI_ROOT}"
+# that will collide with this project's final output! so instead
+TEMPLATE_IOAPI_FN='aqmeii_ag_soil_yearly.nc'
+export TEMPLATE_IOAPI_FP="${WORK_DIR}/${TEMPLATE_IOAPI_FN}"
+export TEMPLATE_IOAPI_DATAVAR_NAME="${AQMEIINA_DV_NAME}"
+export TEMPLATE_IOAPI_TFLAG_NAME="${AQMEIINA_TFLAG_NAME}"
+export TEMPLATE_IOAPI_BAND='1' # since position of TSTEP in N2O(TSTEP, LAY, ROW, COL)? seems to work
+
 ### map scale factors (MSFs)
 # load MSFs from some GRIDCRO2D_<date/> which is pretty large, so gzip it
 
 
 ## global metadata for both {template, "real" monthly/hourly data container} files
 export CLM_CN_TARGET_FILE_ATTR_HISTORY="see ${PROJECT_WEBSITE}"
+export CLM_CN_TARGET_FILE_ATTR_FILEDESC='monthly N2O emissions from natural soils over AQMEII-NA in mol/s, IOAPI-zed'
+# IOAPI metadata: see http://www.baronams.com/products/ioapi/INCLUDE.html#fdesc
+export IOAPI_TIMESTEP_LENGTH='10000' # :TSTEP
+export IOAPI_START_OF_TIMESTEPS='0'  # :STIME
 
 ## metadata for datavar in both {template, "real" monthly/hourly data container} files
 export CLM_CN_TARGET_DATAVAR_NAME="${CLM_CN_REGRID_DATAVAR_NAME}"
   for CMD in \
     'get_regrid_utils' \
     'get_bash_utils' \
+    'get_filepath_funcs' \
     'get_IOAPI_funcs' \
     'get_lat_lon_funcs' \
     'get_stats_funcs' \
 } # end function get_bash_utils
 
 # isa regrid_utils from https://bitbucket.org/tlroche/regrid_utils
+# To override, copy/mod to ${FILEPATH_FUNCS_FP} before running this script.
+function get_filepath_funcs {
+  if [[ -z "${FILEPATH_FUNCS_FP}" ]] ; then
+    echo -e "${THIS_FN}: ERROR: FILEPATH_FUNCS_FP not defined"
+    exit 7
+  fi
+  if [[ ! -r "${FILEPATH_FUNCS_FP}" ]] ; then
+    # copy from downloaded regrid_utils
+    for CMD in \
+      "cp ${REGRID_UTILS_DIR}/${FILEPATH_FUNCS_FN} ${FILEPATH_FUNCS_FP}" \
+    ; do
+      echo -e "$ ${CMD}"
+      eval "${CMD}"
+    done
+  fi
+  if [[ ! -r "${FILEPATH_FUNCS_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: FILEPATH_FUNCS_FP=='${FILEPATH_FUNCS_FP}' not readable"
+    exit 8
+  fi
+} # end function get_filepath_funcs
+
+# isa regrid_utils from https://bitbucket.org/tlroche/regrid_utils
 # To override, copy/mod to ${IOAPI_FUNCS_FP} before running this script.
 function get_IOAPI_funcs {
   if [[ -z "${IOAPI_FUNCS_FP}" ]] ; then
   fi
 } # end function get_conserv
 
+# TODO: remove/rename 'get_template_input' \
 function setup_resources {
   for CMD in \
     'get_raw_input' \
     'get_template_input' \
+    'get_template_IOAPI' \
     'get_regrid' \
     'get_MSFs' \
     'get_target_template' \
   fi
 } # end function get_template_input
 
+### get the "template" file used for IOAPI-writing
+function get_template_IOAPI {
+  if [[ -z "${TEMPLATE_IOAPI_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: TEMPLATE_IOAPI_FP not defined"
+    exit 8
+  fi
+  if [[ ! -r "${TEMPLATE_IOAPI_FP}" ]] ; then
+    if [[ ! -r "${TEMPLATE_IOAPI_GZ_FP}" ]] ; then
+      for CMD in \
+        "${WGET_TO_FILE} ${TEMPLATE_IOAPI_GZ_FP} ${TEMPLATE_IOAPI_GZ_URI}" \
+      ; do
+        echo -e "$ ${CMD}"
+        eval "${CMD}"
+      done
+    fi
+    if [[ -r "${TEMPLATE_IOAPI_GZ_FP}" ]] ; then
+      # JIC ${TEMPLATE_IOAPI_FP} != ${TEMPLATE_IOAPI_GZ_FP} - .gz
+      for CMD in \
+        "gunzip -c ${TEMPLATE_IOAPI_GZ_FP} > ${TEMPLATE_IOAPI_FP}" \
+      ; do
+        echo -e "$ ${CMD}"
+        eval "${CMD}"
+      done
+    fi
+  fi
+  if [[ ! -r "${TEMPLATE_IOAPI_FP}" ]] ; then
+    echo -e "${THIS_FN}::${FUNCNAME[0]}: ERROR: cannot read TEMPLATE_IOAPI_FP=='${TEMPLATE_IOAPI_FP}'"
+    exit 9
+  fi
+} # end function get_template_IOAPI
+
 ### get the regridded output from the R script? no:
 ### only for testing. TODO: provide commandline-switch control for this
 function get_regrid {

check_conservation.ncl

     ; count NAs in raw-input month
     in_nNA(i_month) = num(ismissing(in_data))
   
-; start debugging
-print(str_get_nl() + "in_sum("+i_month+")=="+in_sum(i_month))
-print("in_nNA("+i_month+")=="+in_nNA(i_month))
-;   end debugging
+; ; start debugging
+; print(str_get_nl() + "in_sum("+i_month+")=="+in_sum(i_month))
+; print("in_nNA("+i_month+")=="+in_nNA(i_month))
+; ;   end debugging
 
     ;;; process final outputs
     ;; open monthly output file, noting NCL does not support *.ncf :-(
 
 ;    end do ; iterate hours
   
-; start debugging
-print("out_sum("+i_month+")=="+out_sum(i_month))
-print("out_nNA("+i_month+")=="+out_nNA(i_month))
-;   end debugging
+; ; start debugging
+; print("out_sum("+i_month+")=="+out_sum(i_month))
+; print("out_nNA("+i_month+")=="+out_nNA(i_month))
+; ;   end debugging
       
     ;;; cleanup
     ;; "close" monthly file
 ;!/usr/bin/env ncl ; requires version >= 5.1.1
+;----------------------------------------------------------------------
 ; Convert regridded CLM-CN data to the temporality required for CMAQ.
+; Assumes a bash/linux environment with `cp`, `mv`, `rm`, `ncdump`
 
 ;----------------------------------------------------------------------
 ; libraries
 ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
 ;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; for 
 
-;----------------------------------------------------------------------
-; functions
-;----------------------------------------------------------------------
-
 ; can't call `getenv` before `begin`?
-; can't even load from a string variable?
+load "$FILEPATH_FUNCS_FP"
+load "$IOAPI_FUNCS_FP"
 load "$SUMMARIZE_FUNCS_FP"
 load "$TIME_FUNCS_FP"
 
 ;----------------------------------------------------------------------
+; functions
+;----------------------------------------------------------------------
+
+;----------------------------------------------------------------------
 ; code
 ;----------------------------------------------------------------------
 
   regrid_coordvar_time_name = getenv("CLM_CN_REGRID_TIME_COORDVAR_NAME")
 
   ;; create data container files from extents/template file
-  template_fp = getenv("TEMPLATE_INPUT_FP")
-  template_datavar_name = getenv("TEMPLATE_DATAVAR_NAME")
+  template_fp = getenv("TEMPLATE_IOAPI_FP")
+;  template_datavar_name = getenv("TEMPLATE_DATAVAR_NAME") ; TODO: rename in driver
+  template_datavar_name = getenv("TEMPLATE_IOAPI_DATAVAR_NAME")
+  template_TFLAG_name = getenv("TEMPLATE_IOAPI_TFLAG_NAME")
   ; note template has no coordvars
 
   ; driver must
   output_fp_template_NCL = getenv("CLM_CN_TARGET_FP_TEMPLATE_NCL")
   ; replace output_fp_template_str with <year/><month/> to create "real" filepaths
   output_fp_template_str = getenv("CLM_CN_TARGET_FN_TEMPLATE_STR")
-  out_history_attr = getenv("CLM_CN_TARGET_FILE_ATTR_HISTORY")
+  output_attr_filedesc = getenv("CLM_CN_TARGET_FILE_ATTR_FILEDESC")
+  output_attr_history = getenv("CLM_CN_TARGET_FILE_ATTR_HISTORY")
+  output_attr_stime = getenv("IOAPI_START_OF_TIMESTEPS")
+  output_attr_tstep = getenv("IOAPI_TIMESTEP_LENGTH")
   ; datavar metadata
-  out_datavar_name = getenv("CLM_CN_TARGET_DATAVAR_NAME")
-  out_datavar_attr_long_name = getenv("CLM_CN_TARGET_DATAVAR_ATTR_LONG_NAME")
-  out_datavar_attr_units = getenv("CLM_CN_TARGET_DATAVAR_ATTR_UNITS")
-  out_datavar_attr_var_desc = getenv("CLM_CN_TARGET_DATAVAR_ATTR_VAR_DESC")
-  ; coordvars
-  out_coordvar_x_name = getenv("CLM_CN_TARGET_X_COORDVAR_NAME")
-  out_coordvar_y_name = getenv("CLM_CN_TARGET_Y_COORDVAR_NAME")
-  out_coordvar_layer_name = getenv("CLM_CN_TARGET_LAYER_COORDVAR_NAME")
-  out_coordvar_time_name = getenv("CLM_CN_TARGET_TIME_COORDVAR_NAME")
+  output_dv_name = getenv("CLM_CN_TARGET_DATAVAR_NAME")
+;   out_datavar_attr_long_name = getenv("CLM_CN_TARGET_DATAVAR_ATTR_LONG_NAME")
+;   out_datavar_attr_units = getenv("CLM_CN_TARGET_DATAVAR_ATTR_UNITS")
+;   out_datavar_attr_var_desc = getenv("CLM_CN_TARGET_DATAVAR_ATTR_VAR_DESC")
+;   ; coordvars
+;   out_coordvar_x_name = getenv("CLM_CN_TARGET_X_COORDVAR_NAME")
+;   out_coordvar_y_name = getenv("CLM_CN_TARGET_Y_COORDVAR_NAME")
+;   out_coordvar_layer_name = getenv("CLM_CN_TARGET_LAYER_COORDVAR_NAME")
+;   out_coordvar_time_name = getenv("CLM_CN_TARGET_TIME_COORDVAR_NAME")
 
 ;----------------------------------------------------------------------
 ; payload
 ; system("cp ~/code/regridding/CLM_CN_global_to_AQMEII-NA/2008PTONCLMCNN2O_reunit_regrid.nc ./")
 
 ;----------------------------------------------------------------------
-; create template output
+; write output as monthly, IOAPI
 ;----------------------------------------------------------------------
 
-  ;;; setup template file for monthly output
-  out_fh = addfile(output_fp_template_NCL, "c")
-
-  ;; global attributes
-  ; delete all global attributes? nope, because we cannot delete "file variables":
-  ; http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclVariables.shtml (bold in original)
-  ; > dimensions, coordinate variables and variables cannot be deleted from a file.
-  ; add global attributes={HISTORY, YEAR} now (will add attr=MONTH later)
-  out_fh_att = True ; assign file/global attributes
-  out_fh_att@HISTORY = out_history_attr
-  out_fh_att@YEAR = model_year
-; do once for monthlies
-;  fileattdef(out_fh, out_fh_att) ; copy file attributes   
-
-  ;;; get regridded (and reunit-ed) data
+  ;;; TODO: CACHE! repeating lots of static reads/writes 12 times
+  ;;; setup annual/stable inputs
   regrid_fh = addfile(regrid_fp, "r")
   regrid_dv = regrid_fh->$regrid_datavar_name$
-  regrid_coordvar_x = regrid_dv&$regrid_coordvar_x_name$
-  regrid_coordvar_y = regrid_dv&$regrid_coordvar_y_name$
-  regrid_coordvar_time = regrid_dv&$regrid_coordvar_time_name$
-
-  ;;; get template for output schema
   template_fh = addfile(template_fp, "r")
   template_dv = template_fh->$template_datavar_name$
-;  template has no coordvars
-;  template_coordvar_layer = template_dv&$template_coordvar_layer_name$
-
-  ;;; create datavar and attributes
-  ; rename datavar? no.
-  ; copy/mod datavar? theoretically, but I can't rename or change attributes *on file variables*.
-  ; so instead,
-  ; create a new "datavar" with sizeof(the template datavar) but (eventually) data from the regrid datavar
-  temp_var = new(dimsizes(template_dv), float, getFillValue(regrid_dv))
-  ; note "getFillValue [is] only intended to be used within the new statement."
-
-; ; start debugging---------------------------------------------------
-; print("printVarSummary(temp_var)==")
-; printVarSummary(temp_var)
-; ;   end debugging---------------------------------------------------
-
-  ; coordinate vars
-  temp_var!0 = out_coordvar_time_name ; C-style indexing!
-  temp_var&$out_coordvar_time_name$ = ispan(1,25,1)
-
-  temp_var!1 = out_coordvar_layer_name
-  ; layer "dimensions" can be ignored? template has no coordvars
-;  temp_var&$out_coordvar_layer_name$ = template_coordvar_layer
-
-  temp_var!2 = out_coordvar_y_name
-  temp_var&$out_coordvar_y_name$ = regrid_coordvar_y
-
-  temp_var!3 = out_coordvar_x_name
-  temp_var&$out_coordvar_x_name$ = regrid_coordvar_x
-
-; ; start debugging---------------------------------------------------
-; print("printVarSummary(temp_var)==")
-; printVarSummary(temp_var)
-; ;   end debugging---------------------------------------------------
-
-  ; metadata
-  temp_var@units = out_datavar_attr_units
-  temp_var@long_name = out_datavar_attr_long_name
-  temp_var@var_desc = out_datavar_attr_var_desc
-
-; ; start debugging---------------------------------------------------
-; print("printVarSummary(temp_var)==")
-; printVarSummary(temp_var)
-; ;   end debugging---------------------------------------------------
-
-  ; write the "new datavar" to the output template file
-  out_fh->$out_datavar_name$ = temp_var
-  delete(temp_var) ; since I reuse below
-  out_dv = out_fh->$out_datavar_name$ ; TODO: can I do this first?
-
-; ; start debug---------------------------------------------------------
-; printVarSummary(out_fh)
-; ; Variable: out_fh
-; ; Type: file
-; ; File path:  ./emis_mole_N2O_######_12US1_cmaq_cb05_soa_2008ab_08c.nc
-; ; Number of global attributes:         0
-; ; Number of dimensions:        4
-; ; Number of variables:         4
-
-; printFileVarSummary(out_fh, out_datavar_name)
-; ; Variable: N2O
-; ; Type: float
-; ; Total Size: 13724100 bytes
-; ;             3431025 values
-; ; Number of Dimensions: 4
-; ; Dimensions and sizes:	[TSTEP | 25] x [LAY | 1] x [ROW | 299] x [COL | 459]
-; ; Coordinates: 
-; ;             TSTEP: [1..25]
-; ;             LAY: not a coordinate variable
-; ;             ROW: [1854000..-1722000]
-; ;             COL: [-2550000..2946000]
-; ; Number of Attributes: 4
-; ;   var_desc :	Model species N2O                                                               
-; ;   long_name :	N2O             
-; ;   units :	moles/s         
-; ;   _FillValue :	-3.4e+38
-; ;   end debug---------------------------------------------------------
-
-;----------------------------------------------------------------------
-; create temporalized output
-;----------------------------------------------------------------------
-
-  ;;; retemporalize cleaned (and previously reunit-ed) regrid file: see
-  ;;; Simulation-of-N2O-Production-and-Transport-in-the-US-Cornbelt-Compared-to-Tower-Measurements#wiki-input-processing-CLM-CNv3.5
-  ;;; section="convert concentrations")
-
+  template_TFLAG = template_fh->$template_TFLAG_name$ ; "fake" IOAPI metadatavar
+  
   ;;; write output to monthly files
+;   i_month=0
   do i_month=0 , 11 ; remember: NCL indexing is {C-style,0-based}
     ;; copy template file to monthly file
     n_month = i_month + 1
-    output_fp_monthly_NCL = \
-      str_sub_str(output_fp_template_NCL, output_fp_template_str, \
-        sprinti("%4i", model_year)+sprinti("%0.2i", n_month))
-    output_fp_monthly_CMAQ = \
-      str_sub_str(output_fp_template_CMAQ, output_fp_template_str, \
-        sprinti("%4i", model_year)+sprinti("%0.2i", n_month))
-    output_fp_monthly_NCL_fh = addfile(output_fp_monthly_NCL, "c")
-
-    ; start progress------------------------------------------------------
-    print(str_get_nl() + \
-      this_fn+": processing output from month="+sprinti("%0.2i", n_month)+" to " + \
-      systemfunc("basename "+output_fp_monthly_CMAQ))
-    ;   end progress------------------------------------------------------
-
-    ;; write global attributes to monthly file
-    output_fp_monthly_fh_att = out_fh_att
-    ; TODO: non-system-dependent date
-    output_fp_monthly_fh_att@creation_date = systemfunc("date")
-    output_fp_monthly_fh_att@MONTH = n_month
-    fileattdef( \
-      output_fp_monthly_NCL_fh, \
-      output_fp_monthly_fh_att) ; copy file attributes   
-
-    ;; write datavar definition to monthly file
-;     copy_VarMeta( \
-;       out_fh->$out_datavar_name$, \
-;       output_fp_monthly_NCL_fh->$out_datavar_name$)
-    ; gotta use the temporary datavar trick to copy (again)
-    temp_var = new(dimsizes(out_dv), float, getFillValue(out_dv))
-    copy_VarMeta(out_dv, temp_var)
-    ; ... and then write the "new datavar" to the file, and ...
-; does not exist yet!
-;    output_fp_monthly_NCL_dv = output_fp_monthly_NCL_fh->$out_datavar_name$
-    output_fp_monthly_NCL_fh->$out_datavar_name$ = temp_var
-
-; FAIL! DON'T DO THIS! output_fp_monthly_NCL_dv is separate from the file! writes to it don't write the file!
-;     output_fp_monthly_NCL_dv = output_fp_monthly_NCL_fh->$out_datavar_name$
-
-    ;; ... then write ONLY datavar data (not metadata!) to monthly file: see
-    ;; http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclVariables.shtml#ValueOnly
-    ; ...(TSTEP | hour, LAY | 1, ROW | :, COL | :)
-    monthly_arr = regrid_dv(i_month,:,:)
-
-; ; start debugging
-; print(str_get_nl() + "from regrid_dv:")
-; print("monthly_sum("+i_month+")=="+sum(monthly_arr))
-; print("monthly_nNA("+i_month+")=="+num(ismissing(monthly_arr)))
-; ;   end debugging
-
-    do i_hour = 0 , 24 ; yes, CMAQ days have 25 hours, last is duplicated
-      ;;; temporality: every hour in a month, emit the same molar-mass rate (mol/s)
-;       output_fp_monthly_NCL_dv(i_hour,0,:,:) = monthly_arr ; FAIL! DON'T DO THIS!
-      output_fp_monthly_NCL_fh->$out_datavar_name$(i_hour,0,:,:) = monthly_arr
-    end do
-
-    ;;; rename the monthly file. TODO: make non-system-dependent, or bug NCL to make this unnecessary
-
-; ; start debugging
-; ; why is check_conservation getting only missing values from files?
-; ; monthly_arr_2 = output_fp_monthly_NCL_dv(0,0,:,:) ; any hour should do
-; monthly_arr_2 = output_fp_monthly_NCL_fh->$out_datavar_name$(0,0,:,:) ; any hour should do
-; print("from PRE-closure " + systemfunc("basename "+output_fp_monthly_NCL) +"::"+out_datavar_name+":")
-; print("monthly_sum("+i_month+")=="+sum(monthly_arr_2))
-; print("monthly_nNA("+i_month+")=="+num(ismissing(monthly_arr_2)))
-; ;   end debugging
-
-;; close the file handle: see http://www.ncl.ucar.edu/Support/talk_archives/2012/2196.html
-;     delete(output_fp_monthly_NCL_dv) ; no fix for bad write--see above
-    delete(output_fp_monthly_NCL_fh)
-
-; ; start debugging
-; ; why is check_conservation getting only missing values from files?
-; monthly_fh_3 = addfile(output_fp_monthly_NCL, "r")
-; monthly_dv_3 = monthly_fh_3->$out_datavar_name$
-; monthly_arr_3 = monthly_dv_3(0,0,:,:) ; any hour should do
-; print("from POST-closure " + systemfunc("basename "+output_fp_monthly_NCL) +"::"+out_datavar_name+":")
-; print("monthly_sum("+i_month+")=="+sum(monthly_arr_3))
-; print("monthly_nNA("+i_month+")=="+num(ismissing(monthly_arr_3)))
-; delete(monthly_fh_3)
-; ;   end debugging
-
-    ;; then make file have the suffix we need :-(
-    cmd = "mv "+output_fp_monthly_NCL+" "+output_fp_monthly_CMAQ
-;    print(this_fn+": about to do: '"+cmd+"'") ; debugging
+
+    ;;; see NCL file-writing procedure outline @ http://www.ncl.ucar.edu/Applications/method_2.shtml
+
+    ;;; 1. Setup monthly filepaths, filehandles, datavars
+    ;; inputs setup above
+    output_fp_monthly_NCL = get_year_and_month_fp( \
+      output_fp_template_NCL, output_fp_template_str, model_year, n_month)
+    output_fp_monthly_CMAQ = get_year_and_month_fp( \
+      output_fp_template_CMAQ, output_fp_template_str, model_year, n_month)
+
+    ;; don't overwrite existing output
+    if (isfilepresent(output_fp_monthly_CMAQ)) then
+      print(this_fn+": file="+output_fp_monthly_CMAQ+" exists, skipping")
+    else
+
+      ;; do overwrite existing intermediate
+      if (isfilepresent(output_fp_monthly_NCL)) then
+        cmd = "rm "+output_fp_monthly_NCL ; TODO: less system-dependently?
+        print(this_fn+": about to do `"+cmd+"`")
+        system(cmd)
+      end if
+
+      output_fh = addfile(output_fp_monthly_NCL, "c")
+      ; below will write output datavar, CMAQ file
+
+      ; start progress------------------------------------------------------
+      print(str_get_nl() + \
+        this_fn+": processing output from month="+sprinti("%0.2i", n_month)+" to " + \
+        systemfunc("basename "+output_fp_monthly_CMAQ))
+      ;   end progress------------------------------------------------------
+
+      ;; 1.1. (optional) declare output file definition mode
+      setfileoption(output_fh, "DefineMode", True)
+
+      ;;; 2. Define global/file attributes.
+
+      ;; 2.1. Copy global/file attributes from template: see Example 3 in http://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml
+      ;; will we need to do this separately for each monthly?
+      template_fa_names = getvaratts(template_fh)
+      if (.not. all(ismissing(template_fa_names))) then
+        do i=0, dimsizes(template_fa_names)-1
+          fa_name = template_fa_names(i)
+          output_fh@$fa_name$ = template_fh@$fa_name$
+        end do
+      end if
+
+      ;; 2.2. Overwrite non-template global/file attributes (if any).
+      output_fh@FILEDESC = output_attr_filedesc
+      output_fh@HISTORY = output_attr_history
+      output_fh@creation_date = systemfunc("date")
+      ; TODO: set CDATE, CTIME, WDATE, WTIME
+      output_fh@YEAR = model_year
+      output_fh@MONTH = n_month
+      ; TODO: set SDATE, STIME, TSTEP
+      output_fh@SDATE = stringtoint(get_monthly_SDATE(model_year, n_month))
+      output_fh@STIME = stringtoint(output_attr_stime)
+      output_fh@TSTEP = stringtoint(output_attr_tstep)
+
+      ;;; 3. Define {dimensions, coordinate variables} and their dimensionality
+
+      ;;; 3.1. Copy IOAPI metadatavar=TFLAG dimensions and attributes from template to output.
+      template_TFLAG_type = typeof(template_TFLAG)     ; type as string
+      template_TFLAG_mv = getFillValue(template_TFLAG) ; missing value
+
+      ;;; Major aside: IOAPI does not apparently typically set attr={_FillValue, missing_value} on datavar=TFLAG, which makes template_TFLAG_mv==NA
+      ;;; Below we then propagate that to output_TFLAG, and NCL is good with that.
+      ;;; But when we try to use this with R::ncdf4, it causes problems. So set template_TFLAG_mv, and
+      ;;; TODO: raise issue on netcdfgroup@unidata.ucar.edu: NCL OK with putting datavar with missing _FillValue, ncdf4 is not
+
+      if (ismissing(template_TFLAG_mv)) then
+        template_TFLAG_mv = default_fillvalue(template_TFLAG_type)
+      end if
+      template_TFLAG_dims_sizes = dimsizes(template_TFLAG)
+      template_TFLAG_dims_n = dimsizes(template_TFLAG_dims_sizes) ; == |dims|
+      template_TFLAG_dims_names = getvardims(template_TFLAG) ; TSTEP, VAR, DATE-TIME
+      vars_n = template_TFLAG_dims_sizes(1)
+      datetimes_n = template_TFLAG_dims_sizes(2)
+
+      ;;; 3.2. Copy "real" datavar dimensions from template to output.
+
+      output_dv_type = typeof(template_dv)     ; type as string
+      output_dv_mv = getFillValue(template_dv) ; missing value
+      if (ismissing(output_dv_mv)) then
+        output_dv_mv = default_fillvalue(output_dv_type)
+      end if
+      output_dv_dims_names = getvardims(template_dv)
+      output_dv_dims_sizes = dimsizes(template_dv)
+      tsteps_n = output_dv_dims_sizes(0)
+      layers_n = output_dv_dims_sizes(1)
+      rows_n = output_dv_dims_sizes(2)
+      cols_n = output_dv_dims_sizes(3)
+
+      ;;; 3.3. Copy file dimensions from template to output.
+      output_fh_dims_names = getvardims(template_fh)
+      output_fh_dims_n = dimsizes(output_fh_dims_names)
+      output_fh_dims_sizes = new(output_fh_dims_n, integer, "No_FillValue") ; dimensions must never have NA values
+
+      do i_dim = 0 , output_fh_dims_n-1
+        i_dim_name = output_fh_dims_names(i_dim)
+        if      (i_dim_name .eq. "COL") then
+          output_fh_dims_sizes(i_dim) = cols_n
+        else if (i_dim_name .eq. "DATE-TIME") then
+          output_fh_dims_sizes(i_dim) = datetimes_n
+        else if (i_dim_name .eq. "LAY") then
+          output_fh_dims_sizes(i_dim) = layers_n
+        else if (i_dim_name .eq. "ROW") then
+          output_fh_dims_sizes(i_dim) = rows_n
+        else if (i_dim_name .eq. "TSTEP") then
+          output_fh_dims_sizes(i_dim) = tsteps_n
+        else if (i_dim_name .eq. "VAR") then
+          output_fh_dims_sizes(i_dim) = vars_n
+        end if ; why so many?
+        end if  
+        end if  
+        end if  
+        end if  
+        end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+      end do
+
+      ;; Should any dimensions be unlimited? Here, no.
+      output_fh_dims_unlim = new(output_fh_dims_n, logical, "No_FillValue")
+      do i=0 , output_fh_dims_n-1
+        output_fh_dims_unlim(i) = False
+      end do
+
+      ; at last! define all our dimensions
+      filedimdef(output_fh,         \
+        (/ output_fh_dims_names /), \
+        (/ output_fh_dims_sizes /), \
+        (/ output_fh_dims_unlim /)  \
+      )
+
+      ;;; 4. Create output datavars.
+
+      ;;; 4.1 TFLAG dimensions and attributes
+
+      output_TFLAG =  new(template_TFLAG_dims_sizes, template_TFLAG_type, template_TFLAG_mv)
+      ;; ... copy its dimensions and attributes from template
+      copy_VarMeta(template_TFLAG, output_TFLAG)
+      ;; define variable on file
+      filevardef(output_fh, template_TFLAG_name, template_TFLAG_type, template_TFLAG_dims_names)
+      ;; copy attributes to datavar
+      filevarattdef(output_fh, template_TFLAG_name, output_TFLAG)
+
+      ;;; 4.2 "real" datavar dimensions and attributes
+      output_dv =  new(output_dv_dims_sizes, output_dv_type, output_dv_mv)
+      ;; ... copy its dimensions and attributes from template
+      copy_VarMeta(template_dv, output_dv) ; our dimensions are not a "leftmost subset"
+      ;; define variable on file
+      filevardef(output_fh, output_dv_name, output_dv_type, output_dv_dims_names)
+      ;; copy attributes to datavar
+      filevarattdef(output_fh, output_dv_name, output_dv)
+
+      ;; 4.3. (optional) Exit file definition mode
+      setfileoption(output_fh, "DefineMode", False)
+
+      ;;; 5. Write data to output variable(s)
+      output_fh->$template_TFLAG_name$ = (/ output_TFLAG /)
+
+      ;; ...(TSTEP | hour, LAY | 1, ROW | :, COL | :)
+      do i_hour = 0 , 24 ; yes, CMAQ days have 25 hours, last is duplicated
+        ; since all our hours are the same, write all 25 hours
+;         output_fh->$output_dv_name$(i_hour,0,:,:) = (/ regrid_dv(i_month,:,:) /)
+        ; rows are reversed in VERDI! see equivalent GEIA problem
+        output_fh->$output_dv_name$(i_hour,0,:,:) = (/ regrid_dv(i_month,::-1,:) /)
+      end do
+
+      ;;; 6. Close output filehandle to flush file: see http://www.ncl.ucar.edu/Support/talk_archives/2012/2196.html
+      delete(output_fh)
+
+    end if ; isfilepresent(output_fp_monthly_CMAQ)
+
+; start debug-----------------------------------------------------------
+
+    if ((.not. (isfilepresent(output_fp_monthly_CMAQ))) .and. \
+        (.not. (isfilepresent(output_fp_monthly_NCL)))) then
+      print(this_fn+": ERROR: cannot find either monthly_CMAQ or monthly_NCL output")
+      return ; TODO: abend  
+    else if ((isfilepresent(output_fp_monthly_CMAQ)) .and. \
+        (.not. (isfilepresent(output_fp_monthly_NCL)))) then
+      ; symlink it? nope, VERDI is too dumb for that :-(
+      cmd = "cp "+output_fp_monthly_CMAQ+" "+output_fp_monthly_NCL
+      print(this_fn+": about to do `"+cmd+"`")
+      system(cmd)
+    end if ; why so many?
+    end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+
+    ; ASSERT: isfilepresent(output_fp_monthly_NCL)
+    output_fh = addfile(output_fp_monthly_NCL, "r") ; handle
+
+    ; look at the schema
+;     printVarSummary(output_fh)
+;     printVarSummary(output_fh->$output_dv_name$)
+    cmd = "ncdump -h "+output_fp_monthly_NCL
+    print(this_fn +": about to do: `" +cmd+"`")
     system(cmd)
 
-; ; start debugging
-; ; why is check_conservation getting only missing values from files?
-; monthly_fh_4 = addfile(output_fp_monthly_CMAQ, "r")
-; monthly_dv_4 = monthly_fh_4->$out_datavar_name$
-; monthly_arr_4 = monthly_dv_4(0,0,:,:) ; any hour should do
-; print("from " + systemfunc("basename "+output_fp_monthly_CMAQ) +"::"+out_datavar_name+":")
-; print("monthly_sum("+i_month+")=="+sum(monthly_arr_4))
-; print("monthly_nNA("+i_month+")=="+num(ismissing(monthly_arr_4)))
-; delete(monthly_fh_4)
-; ;   end debugging
-
-    ;;; cleanup
-    if (isfilepresent(output_fp_monthly_NCL)) then
-      system("rm "+output_fp_monthly_NCL)
-    end if
-    ; gotta delete(...) one at a time? no:
-    ; > NclBuildArray: can not combine character or string types with numeric types
-    ; gotta make a list
-    delete( [/ \
-      cmd, \
-      output_fp_monthly_CMAQ, \
-      output_fp_monthly_NCL \
-    /] )
-
-  end do ; iterate months
-  
-  ;;; cleanup/teardown
+    ; look at the data
+    print(str_get_nl() + this_fn +": about to do: '" +\
+      "summarize_annual_4d(output_fh->" +\
+      output_dv_name + ")'")
+    summarize_annual_4d(output_fh->$output_dv_name$)
+
+;     ;;; View output in VERDI. Note:
+;     ;;; * 'False'==foreground process.
+;     ;;; * Backgrounding in console causes something to vomit on the ampersand ending the generated VERDI call:
+;     ;;;   it appends ' | more', which fails (?!?)
+;     ;;;   If so, just copy everything up to the ampersand, and wrap in `system("")
+;     show_VERDI_fast_tile_plot(output_fp_monthly_NCL, output_dv_name, False)
+
+;   end debug-----------------------------------------------------------
+
+; ----------------------------------------------------------------------
+; cleanup
+; ----------------------------------------------------------------------
+
+    if ((.not. (isfilepresent(output_fp_monthly_CMAQ))) .and. \
+        (.not. (isfilepresent(output_fp_monthly_NCL)))) then
+      print(this_fn+": ERROR: cannot find either monthly_CMAQ or monthly_NCL output")
+      return ; TODO: abend  
+    else if ((.not. (isfilepresent(output_fp_monthly_CMAQ))) .and. \
+             (isfilepresent(output_fp_monthly_NCL))) then
+      ;; make file have the suffix we need :-( TODO: less system-dependently?
+      cmd = "mv "+output_fp_monthly_NCL+" "+output_fp_monthly_CMAQ
+;       print(this_fn+": about to do: '"+cmd+"'") ; debugging
+      system(cmd)
+    else if ((isfilepresent(output_fp_monthly_CMAQ)) .and. \
+             (isfilepresent(output_fp_monthly_NCL))) then
+      ;; remove file we don't need :-( TODO: less system-dependently?
+      cmd = "rm "+output_fp_monthly_NCL
+;       print(this_fn+": about to do: '"+cmd+"'") ; debugging
+      system(cmd)
+    end if ; why so many?
+    end if ; see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
+    end if ; isfilepresent(output_fp_monthly_CMAQ)
+    ; ASSERT: ((isfilepresent(output_fp_monthly_CMAQ)) .and.
+    ;          (.not. (isfilepresent(output_fp_monthly_NCL))))
+
+  end do ; i_month=0 , 11
 
 end ; retemp.ncl