Tom Roche avatar Tom Roche committed baa15ec

fix grid alignment, fully IOAPIize hourly output, enable view in VERDI

After fully IOAPI-ize output from make_hourlies.ncl, and viewing in VERDI (which also verifies goodness of IOAPI-ization), I was able to
* verify that make_hourlies.ncl's output was flipped (just like the previous projects')
* fix it

Also:

* add lots debugging and progress code to make_hourlies.ncl
* add lots progress code to check_conservation.ncl (more below)
* removed no-longer-used output code from make_hourlies.ncl
* added a few more generic IOAPI-writing strings (needed for make_hourlies.ncl) to GFED_driver.sh (TODO: refactor to uber-project)

Unfortunately, even running 'env GIT_SSL_NO_VERIFY=true ./GFED_driver.sh' on terrae/EMVL in fresh terminal takes too long, due to check_conservation.ncl
(first month of aggregating hourly takes 18 min wallclock) so I'm committing hopefully/foolishly while it runs.
TODO: run uber_driver.sh:repo overnight

TODO:
* move all these TODOs to https://bitbucket.org/tlroche/gfed-3.1_global_to_aqmeii-na/issues
*** investigate 8% loss of mass from monthly regrids to hourly regrids
* add switch to gate running of check_conservation.ncl::aggregation of hourlies: takes too long!
* add switch to gate VERDI-running
** ... and handle {name, path to} VERDI in regrid_utils/bash_utilities.sh
* check_raw_fractions.ncl: check that
ftp://gfed3:dailyandhourly@zea.ess.uci.edu/GFEDv3.1/Readme.pdf
> Over a day at each grid cell location, the sum of the eight 3-hourly fire fractions should be equal to 1.0.
* check_coord_vars.ncl: should also check monotonicity of coordvars? or is that guaranteed by NetCDF?
* 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-style 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.*(...)
* 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)

 ## regrid constants
 export GFED_REGRID_N_ROW='299'
 export GFED_REGRID_N_COL='459'
+# 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
 
 ### numerical and math constants
 export ZERO_MAX='0.0001' # any value less than this is probably 0
 export GFED_REGRID_HOURLY_EMIS_JAN_01_FP_NCL="${GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL}/${GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING}/0101}"
 
 ## file/global attributes
-export GFED_REGRID_HOURLY_EMIS_HISTORY="see ${PROJECT_WEBSITE}"
-export GFED_REGRID_HOURLY_EMIS_TITLE="hourly N2O emissions due to biomass burning for ${MODEL_YEAR}, derived from GFED emissions and fractions"
+export GFED_REGRID_HOURLY_EMIS_ATTR_HISTORY="see ${PROJECT_WEBSITE}"
+export GFED_REGRID_HOURLY_EMIS_ATTR_TITLE='hourly N2O emissions due to biomass burning derived from GFED emissions and fractions over AQMEII-NA in mol/s, IOAPI-zed'
+export GFED_REGRID_HOURLY_EMIS_ATTR_FILEDESC="${GFED_REGRID_HOURLY_EMIS_ATTR_TITLE}"
 
 ### datavar
 export GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME='N2O'

check_conservation.ncl

 ; process (regridded) hourlies
 ;----------------------------------------------------------------------
 
+  print(str_get_nl() + \
+    this_fn+": about to aggregate hourly emissions--will take awhile!" + str_get_nl() + \
+    "starting aggregation @ " + systemfunc("date") + str_get_nl()) ; TODO: better progress control
+
   ;; declare output
   hourly_arr = new(hourly_arr_dims, "float")
 
-  print(str_get_nl() + \
-    this_fn+": about to aggregate hourly emissions--will take awhile!" + \
-  str_get_nl()) ; TODO: better progress control
-
   ;; loop over files
 ;  do n_month = 1, 2 ; for testing
 ;  do n_month = 12, 12 ; for debugging
   do n_month = 1, 12
+    print(str_get_nl() + \
+      this_fn+": aggregating hourly emissions for month=" + n_month + str_get_nl())
     i_month = n_month - 1 ; array pointer
     ;;; write the hourly emissions for the day
     first_day_this_month = \ ; will need this later
       first_day_of_month(model_year, n_month)
     do i_day = first_day_this_month, last_day_of_month(model_year, n_month)
-      ; TODO: show progress
       ;; i_day is sequential/Julian, so don't use that
       n_day = i_day - first_day_this_month + 1
+      print(str_get_nl() + \
+        this_fn+": aggregating hourly emissions for " + \
+        sprinti("%0.2i", n_month) + sprinti("%0.2i", n_day) + " @ " + \
+        systemfunc("date") + str_get_nl())
+
       hourly_fp = get_daily_fp(\
         hourly_fp_template, hourly_template_string, n_month, n_day)
       hourly_fh = addfile(hourly_fp, "r")
   ; start progress----------------------------------------------------
   ; TODO: better progress control
   print(str_get_nl() + \
-    this_fn+": start summarizing aggregate of hourlies @ " + \
+    this_fn+": about to summarize aggregate of hourlies @ " + \
     systemfunc("date") + str_get_nl() + \
     "will take awhile! (15-90 min on terrae interactive)" + str_get_nl() )
   ;   end progress----------------------------------------------------

make_hourlies.ncl

 
 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 
+load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; for getFillValue
 
 ; can't call `getenv` before `begin`?
 load "$FILEPATH_FUNCS_FP"  ; for get_daily_fp
 load "$IOAPI_FUNCS_FP"     ; IOAPI-writing and testing
-load "$STRING_FUNCS_FP"    ; for right_pad_blanks
+load "$SUMMARIZE_FUNCS_FP"    ; for summarize_array
 load "$TIME_FUNCS_FP"      ; for first_day_of_month, last_day_of_month
 
 ;----------------------------------------------------------------------
 hours_per_day = stringtoint(getenv("HOURS_PER_DAY"))
 seconds_per_hour = stringtofloat(getenv("SECONDS_PER_HOUR"))
 molar_mass_N2O = stringtofloat(getenv("MOLAR_MASS_N2O"))
+;; visualization
+sigdigs = stringtoint(getenv("OUTPUT_SIGNIFICANT_DIGITS"))
 
 ;;; inputs
 
 hour3ly_frac_dv_name = getenv("GFED_REGRID_3HOURLY_FRAC_DATAVAR_NAME")
 
 ;; AQMEII-NA gridcell areas: don't actually need here (?) but am leaving in for future reference
-out_areas_fp = getenv("AQMEIINA_AREAS_FP")
-out_areas_dv_name = getenv("AQMEIINA_AREAS_DATAVAR_NAME")
+; out_areas_fp = getenv("AQMEIINA_AREAS_FP")
+; out_areas_dv_name = getenv("AQMEIINA_AREAS_DATAVAR_NAME")
+
+;; create data container files from extents/template file
+template_fp = getenv("TEMPLATE_INPUT_FP")
+template_dv_name = getenv("TEMPLATE_INPUT_DATAVAR_NAME")
+template_TFLAG_name = getenv("TEMPLATE_INPUT_TFLAG_NAME")
+; note template has no coordvars
 
 ;;; outputs
 
 out_emis_fp_template_cmaq = getenv("GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_CMAQ")
 out_emis_fp_template_ncl = getenv("GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL")
 ; attributes: TODO: port function get_output_file_attributes
-out_attr_file_history = getenv("GFED_REGRID_HOURLY_EMIS_HISTORY")
+out_attr_history = getenv("GFED_REGRID_HOURLY_EMIS_ATTR_HISTORY")
+out_attr_filedesc = getenv("GFED_REGRID_HOURLY_EMIS_ATTR_FILEDESC")
+out_attr_stime = stringtoint(getenv("IOAPI_START_OF_TIMESTEPS"))
+out_attr_tstep = stringtoint(getenv("IOAPI_TIMESTEP_LENGTH"))
 
 ;; datavar
 out_dv_name = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME")
 ; functions
 ;----------------------------------------------------------------------
 
-;;; Create attributes for output data variable.
-;;; Values below copy/mod from FIXME.
-undef("get_output_var_attributes") ; allow redefinition in console
-;  dv [*] [*] [*] [*] : float ; datavar
-function get_output_var_attributes(
-  dv : numeric ; datavar
-)
-begin
-  dv@units = out_dv_attr_units
-  dv@long_name = out_dv_attr_long_name
-  dv@var_desc = out_dv_attr_var_desc
-  return(dv)
-end ; function get_output_var_attributes
-
 ;;; Create file/global attributes for output IOAPI.
 ;;; Values below copied mostly verbatim from FIXME.
 undef("get_output_file_attributes") ; allow redefinition in console
   ; why did IOAPI put a dash in there?
   out_file_att@$var_list_name$ = "N2O             "
 ; out_file_att@FILEDESC = ; loonnngggg string, omitted
-  out_file_att@HISTORY = out_attr_file_history
+  out_file_att@HISTORY = out_attr_history
   return(out_file_att)
 end ; function get_output_file_attributes
 
 ; payload
 ;----------------------------------------------------------------------
 
-  ;;; open regridded emissions and fractions
+  ;;; open regridded emissions, fractions, and template
   ;; monthly emissions
   monthly_emis_fh = addfile(monthly_emis_fp, "r")
   monthly_emis_dv = monthly_emis_fh->$monthly_emis_dv_name$
-;   monthly_emis_arr = (/ monthly_emis_dv /)
-;   monthly_emis_dims = dimsizes(monthly_emis_arr) ; (/ 12, 299, 459 /)
+  ;; reverse rows to avoid problems below
+  monthly_emis_dv = monthly_emis_dv(:,::-1,:)
 
   ;; daily fractions (by day for year)
   daily_frac_fh = addfile(daily_frac_fp, "r")
   daily_frac_dv = daily_frac_fh->$daily_frac_dv_name$
-;  daily_frac_arr = (/ daily_frac_dv /)
+  ;; reverse rows to avoid problems below
+  daily_frac_dv = daily_frac_dv(:,::-1,:)
 
   ;; 3-hourly fractions (by month for year)
   hour3ly_frac_fh = addfile(hour3ly_frac_fp, "r")
   hour3ly_frac_dv = hour3ly_frac_fh->$hour3ly_frac_dv_name$
-;  hour3ly_frac_arr = (/ hour3ly_frac_dv /)
-
-  ;;; setup daily emissions output
-  ;; create output attributes
-  out_file_atts = get_output_file_attributes()
-;   ; start debugging---------------------------------------------------
-;   printVarSummary(out_file_atts)
-;   ;   end debugging---------------------------------------------------
-
-  ;; Create output datavar with dimensions=(TSTEP | 25, LAY | 1, ROW | 299, COL | 459)
-  ;; using 3hourly fractions as a template (since that has very similar schema)
-  out_dv = get_output_datavar(\
-    hour3ly_frac_dv, hour3ly_frac_dv_name, hour3ly_frac_fh)
-
-;   ; start debugging---------------------------------------------------
-;   printVarSummary(out_dv) ; check dimensions
-; ; Variable: out_dv
-; ; 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: [1..1]
-; ;             ROW: [1854000..-1722000]
-; ;             COL: [-2550000..2946000]
-; ; Number Of Attributes: 4
-; ;   _FillValue :  9.96921e+36
-; ;   units :       moles/s         
-; ;   long_name :   N2O             
-; ;   var_desc :    N2O                                                                             
-; ; TODO: why are ROW decreasing? and does it matter?
-;   ;   end debugging---------------------------------------------------
+  ;; reverse rows to avoid problems below
+  hour3ly_frac_dv = hour3ly_frac_dv(:,:,::-1,:)
+
+  ;;; template for IOAPI-writing
+  template_fh = addfile(template_fp, "r")
+  template_fa_names = getvaratts(template_fh)
+  template_dv = template_fh->$template_dv_name$
+  template_TFLAG = template_fh->$template_TFLAG_name$ ; "fake" IOAPI metadatavar
+
+  ;; Copy IOAPI metadatavar=TFLAG dimensions and attributes from template
+  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)
+
+  ;; Copy "real" datavar dimensions from template
+  out_dv_type = typeof(template_dv)     ; type as string
+  out_dv_mv = getFillValue(template_dv) ; missing value
+  if (ismissing(out_dv_mv)) then
+    out_dv_mv = default_fillvalue(out_dv_type)
+  end if
+  out_dv_dims_names = getvardims(template_dv)
+  out_dv_dims_sizes = dimsizes(template_dv)
+  tsteps_n = out_dv_dims_sizes(0)
+  layers_n = out_dv_dims_sizes(1)
+  rows_n = out_dv_dims_sizes(2)
+  cols_n = out_dv_dims_sizes(3)
 
   ;;; create unit-conversion matrix to finish mapping
   ;;; from GFED-supplied 3-hourly output (units==gN2O/m^2/(3-hour) per ftp://zea.ess.uci.edu/GFEDv3.1/Readme.pdf )
 ;   out_conv_mx = (/ out_areas_dv /) / (molar_mass_N2O * seconds_per_hour)
 
 ;----------------------------------------------------------------------
-; main loop(s): month, day, hour
+; main loop(s): write output as hourly
 ;----------------------------------------------------------------------
 
   ;;; for each month in *_monthly_emis_emissions_regrid.nc
-;  do n_month = 1, 2 ; for testing
-;  do n_month = 12, 12 ; for debugging
+;   do n_month = 1, 1 ; for testing
   do n_month = 1, 12
     ; TODO: show progress
     i_month = n_month - 1 ; array pointer
+
     ;;; get the month's emissions
+;     ; start debugging---------------------------------------------------
+;     monthly_emis_arr = monthly_emis_dv(i_month,:,:)
+;     printVarSummary(monthly_emis_arr)
+; ; Variable: monthly_emis_arr
+; ; Type: float
+; ; Total Size: 548964 bytes
+; ;             137241 values
+; ; Number of Dimensions: 2
+; ; Dimensions and sizes:   [ROW | 299] x [COL | 459]
+; ; Coordinates: 
+; ;             ROW: [-1722000..1854000]
+; ;             COL: [-2550000..2946000]
+; ; Number Of Attributes: 9
+; ;   max : <ARRAY of 12 elements>
+; ;   min : <ARRAY of 12 elements>
+; ;   projection_format :   PROJ.4
+; ;   projection :  +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
+; ;   long_name :   N2O emissions
+; ;   missing_value :       -3.4e+38
+; ;   units :       molN20/mo
+; ;   TSTEP :       1
+; ;   _FillValue :  -3.4e+38
+;     ;   end debugging---------------------------------------------------
     monthly_emis_arr = (/ monthly_emis_dv(i_month,:,:) /)
-    ; debugging: check dimensions == (/ 299, 459 /)
-    ;;; get the month's 3hourly fractions
-    hour3ly_frac_arr = (/ hour3ly_frac_dv(i_month,:,:,:) /)
-    ; debugging: check dimensions == (/ 8, 299, 459 /)
-    ;;; create "hourly fractions" array for the month
-    hourly_frac_arr = get_hourly_frac_arr(hour3ly_frac_arr)
 
-;     ; start debugging-------------------------------------------------
-;     printVarSummary(hourly_frac_arr) ; check dimensions == (/ 24, 299, 459 /)
-; ; Variable: hourly_frac_arr
+    ;;; get the month's 3hourly fractions
+;     ; start debugging---------------------------------------------------
+;     hour3ly_frac_arr = hour3ly_frac_dv(i_month,:,:,:)
+;     printVarSummary(hour3ly_frac_arr)
+; ; Variable: hour3ly_frac_arr
 ; ; Type: float
-; ; Total Size: 13175136 bytes
-; ;             3293784 values
+; ; Total Size: 4391712 bytes
+; ;             1097928 values
 ; ; Number of Dimensions: 3
-; ; Dimensions and sizes: [24] x [299] x [459]
+; ; Dimensions and sizes:   [LAY | 8] x [ROW | 299] x [COL | 459]
 ; ; Coordinates: 
-; ; Number Of Attributes: 1
-; ;   _FillValue :        9.96921e+36
-;     ;   end debugging-------------------------------------------------
+; ;             LAY: [0..21]
+; ;             ROW: [-1722000..1854000]
+; ;             COL: [-2550000..2946000]
+; ; Number Of Attributes: 4
+; ;   TSTEP :       1
+; ;   _FillValue :  9.96921e+36
+; ;   long_name :   3-hourly fraction of daily N2O emission
+; ;   units :       unitless
+;     ;   end debugging---------------------------------------------------
+    hour3ly_frac_arr = (/ hour3ly_frac_dv(i_month,:,:,:) /)
+
+    ;;; compute "hourly fractions" array for the month
+    hourly_frac_arr = get_hourly_frac_arr(hour3ly_frac_arr)
 
     ;;; for each day in the month in *_daily_frac_fractions_regrid.nc
     ;;; (encoded as sequential int's, i.e., 1 Feb -> 32),
     ;;; write the hourly emissions for the day
     first_day_this_month = \ ; will need this later
       first_day_of_month(model_year, n_month)
-;     last_day_this_month = \  ; debugging
-;       last_day_of_month(model_year, n_month)
+;     do i_day = first_day_this_month, first_day_this_month ; debugging
     do i_day = first_day_this_month, last_day_of_month(model_year, n_month)
 
       ;;; Setup daily emissions file path, or rather paths:
       ;;; to get CMAQ extension, gotta work around {CMAQ odd practice, NCL supported extensions}
-      n_day = i_day - first_day_this_month + 1
+      n_day = i_day - first_day_this_month + 1 ; 1-based day of month 
+;       print(this_fn+": debugging: n_month='"+n_month+"', n_day='"+n_day+"'")
+      sdate = get_SDATE(model_year, n_month, n_day) ; reuse later
+;       print(this_fn+": debugging: SDATE='"+sdate+"'")
       out_fp_cmaq = get_daily_fp(\
         out_emis_fp_template_cmaq, day_template_string, n_month, n_day)
       out_fp_ncl = get_daily_fp(\
         out_emis_fp_template_ncl, day_template_string, n_month, n_day)
+
       ;;; don't remake file if already present
       if      (isfilepresent(out_fp_cmaq)) then
         print(this_fn+": skipping output file='"+out_fp_cmaq+"'")
 
       ;; else, 
       ;; get daily fractions
-      daily_frac_arr = (/ daily_frac_dv((i_day-1),:,:) /) ; C-style indexing!
 ;       ; start debugging-----------------------------------------------
+;       daily_frac_arr = daily_frac_dv((i_day-1),:,:) ; C-style indexing!
 ;       printVarSummary(daily_frac_arr) ; check dimensions
 ; ; Variable: daily_frac_arr
 ; ; Type: float
 ; ; Total Size: 548964 bytes
 ; ;             137241 values
 ; ; Number of Dimensions: 2
-; ; Dimensions and sizes: [299] x [459]
+; ; Dimensions and sizes:   [ROW | 299] x [COL | 459]
 ; ; Coordinates: 
-; ; Number Of Attributes: 1
-; ;   _FillValue :        -3.4e+38
+; ;             ROW: [-1722000..1854000]
+; ;             COL: [-2550000..2946000]
+; ; Number Of Attributes: 9
+; ;   max : <ARRAY of 366 elements>
+; ;   min : <ARRAY of 366 elements>
+; ;   projection_format :   PROJ.4
+; ;   projection :  +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
+; ;   long_name :   daily fraction of monthly N2O emission
+; ;   missing_value :       -3.4e+38
+; ;   units :       unitless
+; ;   TSTEP :       1
+; ;   _FillValue :  -3.4e+38
 ;       ;   end debugging-----------------------------------------------
+      daily_frac_arr = (/ daily_frac_dv((i_day-1),:,:) /) ; C-style indexing!
 
       ;; compute daily emissions
       daily_emis_arr = monthly_emis_arr * daily_frac_arr
-;       ; start debugging-----------------------------------------------
-;       printVarSummary(daily_emis_arr) ; check dimensions
-; ; Variable: daily_emis_arr
-; ; Type: float
-; ; Total Size: 548964 bytes
-; ;             137241 values
-; ; Number of Dimensions: 2
-; ; Dimensions and sizes: [299] x [459]
-; ; Coordinates: 
-; ; Number Of Attributes: 1
-; ;   _FillValue :        -3.4e+38
-;       ;   end debugging-----------------------------------------------
 
+;       ;; for each hour in the day, write hourly emissions to datavar
+;       do n_hour = 1, hours_per_day
+;         ; TODO: show progress
+;         i_hour = n_hour - 1
+;         ;; compute hourly emissions in molN2O/s : note conversions above
+;         hourly_emis_arr = daily_emis_arr * hourly_frac_arr(i_hour,:,:)
+;         ;; write daily emissions datavar (assuming |LAY|==1)
+; ;        out_dv(i_hour,0,:,:) = hourly_emis_arr
+; ; whacks {dimension names, coordinate variables} ROW, COL
+;         out_dv(i_hour,0,:,:) = (/ hourly_emis_arr /)
+;       end do ; n_hour = 1, hours_per_day
+
+;       ;; since CMAQ wants one more TSTEP than actual hours, duplicate last hour
+;       out_dv(hours_per_day,0,:,:) = (/ out_dv(hours_per_day-1,0,:,:) /)
+; ;       ; start debugging-----------------------------------------------
+; ;       printVarSummary(out_dv) ; check dimensions
+; ; ; Variable: out_dv
+; ; ; 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: [1..1]
+; ; ;             ROW: [-1722000..1854000]
+; ; ;             COL: [-2550000..2946000]
+; ; ; Number Of Attributes: 4
+; ; ;   _FillValue :        9.96921e+36
+; ; ;   units :     moles/s         
+; ; ;   long_name : N2O             
+; ; ;   var_desc :  N2O                                                                             
+; ;       ;   end debugging-----------------------------------------------
+
+;----------------------------------------------------------------------
+; write fully-IOAPI-zed output (and verify)
+;----------------------------------------------------------------------
+
+      ;;; see NCL file-writing procedure outline @ http://www.ncl.ucar.edu/Applications/method_2.shtml
+
+      ;;; 1. Setup input and output filepaths, filehandles, datavars
+      ;; mostly setup above
+
+      ;; don't overwrite existing output
+      if (isfilepresent(out_fp_cmaq)) then
+        print(this_fn+": file="+out_fp_cmaq+" exists, skipping")
+        continue ; to next day
+      end if
+      ;; else
+
+      ;; do overwrite existing intermediate
+      if (isfilepresent(out_fp_ncl)) then
+        cmd = "rm "+out_fp_ncl ; TODO: less system-dependently?
+        print(this_fn+": about to do `"+cmd+"`")
+        system(cmd)
+      else
+        print(str_get_nl() + "writing data to out_fp='"+out_fp_ncl+"'")
+      end if
+      out_fh = addfile(out_fp_ncl, "c")
+      ; below will write output datavar, CMAQ file
+
+      ;; 1.1. (optional) declare output file definition mode
+      setfileoption(out_fh, "DefineMode", True)
+
+      ;;; 2. Define global/file attributes.
+;       print(str_get_nl() + "2. Define global/file attributes.") ; debugging
+
+      ;; 2.1. Copy global/file attributes from template: see Example 3 in http://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml
+;       print(str_get_nl() + "2.1. Copy global/file attributes from template: see Example 3 in http://www.ncl.ucar.edu/Document/Functions/Built-in/addfile.shtml") ; debugging
+      if (.not. all(ismissing(template_fa_names))) then
+        do i=0, dimsizes(template_fa_names)-1
+          fa_name = template_fa_names(i)
+          out_fh@$fa_name$ = template_fh@$fa_name$
+        end do
+      end if
+
+      ;; 2.2. Overwrite non-template global/file attributes (if any).
+;       print("2.2. Overwrite non-template global/file attributes (str_get_nl() + if any).") ; debugging
+      out_fh@FILEDESC = out_attr_filedesc
+      out_fh@HISTORY = out_attr_history
+      out_fh@creation_date = systemfunc("date")
+      ; TODO: set CDATE, CTIME, WDATE, WTIME
+      out_fh@YEAR = model_year
+      out_fh@MONTH = n_month
+      out_fh@DAY = n_day
+;       out_fh@SDATE = stringtoint(get_SDATE(model_year, n_month, n_day))
+;       sdate = get_SDATE(model_year, n_month, n_day) ; moved up
+      out_fh@SDATE = stringtoint(sdate)
+      out_fh@STIME = out_attr_stime
+      out_fh@TSTEP = out_attr_tstep
+
+      ;;; 3. Define {dimensions, coordinate variables} and their dimensionality
+;       print(str_get_nl() + "3. Define {dimensions, coordinate variables} and their dimensionality") ; debugging
+
+      ;;; 3.1. Copy IOAPI metadatavar=TFLAG dimensions and attributes from template.
+      ;; do this above, outside of loops!
+
+      ;;; 3.2. Copy "real" datavar dimensions from template.
+      ;; do this above, outside of loops!
+
+      ;;; 3.3. Copy file dimensions from template to output.
+;       print(str_get_nl() + "3.3. Copy file dimensions from template to output.") ; debugging
+      out_fh_dims_names = getvardims(template_fh)
+      out_fh_dims_n = dimsizes(out_fh_dims_names)
+      out_fh_dims_sizes = new(out_fh_dims_n, integer, "No_FillValue") ; dimensions must never have NA values
+
+      do i_dim = 0 , out_fh_dims_n-1
+        i_dim_name = out_fh_dims_names(i_dim)
+        if      (i_dim_name .eq. "COL") then
+          out_fh_dims_sizes(i_dim) = cols_n
+        else if (i_dim_name .eq. "DATE-TIME") then
+          out_fh_dims_sizes(i_dim) = datetimes_n
+        else if (i_dim_name .eq. "LAY") then
+          out_fh_dims_sizes(i_dim) = layers_n
+        else if (i_dim_name .eq. "ROW") then
+          out_fh_dims_sizes(i_dim) = rows_n
+        else if (i_dim_name .eq. "TSTEP") then
+          out_fh_dims_sizes(i_dim) = tsteps_n
+        else if (i_dim_name .eq. "VAR") then
+          out_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.
+      out_fh_dims_unlim = new(out_fh_dims_n, logical, "No_FillValue")
+      do i=0 , out_fh_dims_n-1
+        out_fh_dims_unlim(i) = False
+      end do
+
+      ; at last! define all our dimensions
+      filedimdef(out_fh,         \
+        (/ out_fh_dims_names /), \
+        (/ out_fh_dims_sizes /), \
+        (/ out_fh_dims_unlim /)  \
+      )
+
+      ;;; 4. Create output datavars.
+
+      ;;; 4.1 TFLAG dimensions and attributes
+
+      out_TFLAG =  new(template_TFLAG_dims_sizes, template_TFLAG_type, template_TFLAG_mv)
+      ;; ... copy its dimensions and attributes from template
+      copy_VarMeta(template_TFLAG, out_TFLAG)
+      ;; define variable on file
+      filevardef(out_fh, template_TFLAG_name, template_TFLAG_type, template_TFLAG_dims_names)
+      ;; copy attributes to datavar
+      filevarattdef(out_fh, template_TFLAG_name, out_TFLAG)
+
+      ;;; 4.2 "real" datavar dimensions and attributes
+      out_dv =  new(out_dv_dims_sizes, out_dv_type, out_dv_mv)
+      ;; ... copy its dimensions and attributes from template
+      copy_VarMeta(template_dv, out_dv) ; our dimensions are not a "leftmost subset"
+      ;; define variable on file
+      filevardef(out_fh, out_dv_name, out_dv_type, out_dv_dims_names)
+      ;; copy attributes to datavar
+      filevarattdef(out_fh, out_dv_name, out_dv)
+
+      ;; 4.3. (optional) Exit file definition mode
+      setfileoption(out_fh, "DefineMode", False)
+
+      ;;; 5. Write data to output variable(s)
+      ;; TFLAG
+      out_fh->$template_TFLAG_name$ = (/ out_TFLAG /)
+
+      ;; real data
       ;; for each hour in the day, write hourly emissions to datavar
       do n_hour = 1, hours_per_day
         ; TODO: show progress
         i_hour = n_hour - 1
         ;; compute hourly emissions in molN2O/s : note conversions above
         hourly_emis_arr = daily_emis_arr * hourly_frac_arr(i_hour,:,:)
-          
-;         ; start debugging-----------------------------------------------
-;         printVarSummary(hourly_emis_arr) ; check dimensions
-; ; Variable: hourly_emis_arr
-; ; Type: float
-; ; Total Size: 548964 bytes
-; ;             137241 values
-; ; Number of Dimensions: 2
-; ; Dimensions and sizes:	[299] x [459]
-; ; Coordinates: 
-; ; Number Of Attributes: 1
-; ;   _FillValue :	-3.4e+38
-;         ;   end debugging-----------------------------------------------
-
         ;; write daily emissions datavar (assuming |LAY|==1)
-;        out_dv(i_hour,0,:,:) = hourly_emis_arr
+;        out_fh->$out_dv_name$(i_hour,0,:,:) = hourly_emis_arr
 ; whacks {dimension names, coordinate variables} ROW, COL
-        out_dv(i_hour,0,:,:) = (/ hourly_emis_arr /)
-      end do ; iterating hours
+        out_fh->$out_dv_name$(i_hour,0,:,:) = (/ hourly_emis_arr /)
+      end do ; n_hour = 1, hours_per_day
 
       ;; since CMAQ wants one more TSTEP than actual hours, duplicate last hour
-      out_dv(hours_per_day,0,:,:) = (/ out_dv(hours_per_day-1,0,:,:) /)
+      out_fh->$out_dv_name$(hours_per_day,0,:,:) = \
+        (/ daily_emis_arr * hourly_frac_arr(hours_per_day-1,:,:) /)
 
 ;       ; start debugging-----------------------------------------------
 ;       printVarSummary(out_dv) ; check dimensions
 ; ; Total Size: 13724100 bytes
 ; ;             3431025 values
 ; ; Number of Dimensions: 4
-; ; Dimensions and sizes: [TSTEP | 25] x [LAY | 1] x [ROW | 299] x [COL | 459]
+; ; Dimensions and sizes:   [TSTEP | 25] x [LAY | 1] x [ROW | 299] x [COL | 459]
 ; ; Coordinates: 
-; ;             TSTEP: [1..25]
-; ;             LAY: [1..1]
-; ;             ROW: [1854000..-1722000]
-; ;             COL: [-2550000..2946000]
 ; ; Number Of Attributes: 4
-; ;   _FillValue :        9.96921e+36
-; ;   units :     moles/s         
-; ;   long_name : N2O             
-; ;   var_desc :  N2O                                                                             
+; ;   var_desc :    Model species N2O                                                               
+; ;   long_name :   N2O             
+; ;   _FillValue :  -9.999e+36
+; ;   units :       moles/s         
+;       print(this_fn+": |obs("+sdate+":out_dv)|=="+get_n_obs( (/ out_fh->$out_dv_name$ /) ))
 ;       ;   end debugging-----------------------------------------------
 
-      ;;; Write output daily emissions file to NCL filepath
-      out_fh = addfile(out_fp_ncl, "c") ; create, not write
-      ;; write global attributes: must do this before writing datavar
-      out_file_atts@creation_date = systemfunc("date")
-      fileattdef(out_fh, out_file_atts)
-      ;; write datavar
-      ; TODO: replace with real progress control!
-      print("writing data to out_fp='"+out_fp_ncl+"'") 
-      out_fh->$out_dv_name$ = out_dv
-      ;; replace with desired CMAQ filepath
-      system("mv "+out_fp_ncl+" "+out_fp_cmaq)
-
-    end do ; iterating days
-  end do ; iterating months
+      ;;; 6. Close output filehandle to flush file: see http://www.ncl.ucar.edu/Support/talk_archives/2012/2196.html
+      delete(out_fh)
+      
+; ----------------------------------------------------------------------
+; debugging: visualize data
+; ----------------------------------------------------------------------
+
+      ;;; don't look at every day
+      if (n_day .eq. 1) then
+
+        if ((.not. (isfilepresent(out_fp_cmaq))) .and. \
+            (.not. (isfilepresent(out_fp_ncl)))) then
+          print(this_fn+": ERROR: cannot find either monthly_CMAQ or monthly_NCL output")
+          return ; TODO: abend  
+          else if ((isfilepresent(out_fp_cmaq)) .and. \
+                   (.not. (isfilepresent(out_fp_ncl)))) then
+            ; symlink it? nope, VERDI is too dumb for that :-(
+            cmd = "cp "+out_fp_cmaq+" "+out_fp_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(out_fp_ncl)
+        out_fh = addfile(out_fp_ncl, "r") ; handle
+        out_dv = out_fh->$out_dv_name$
+        out_arr = (/ out_dv /)
+
+        ;;; look at the schema
+;         printVarSummary(out_fh)
+        cmd = "ncdump -h "+out_fp_ncl
+        print(this_fn +": about to do: `" +cmd+"`")
+        system(cmd)
+        printVarSummary(out_dv)
+
+        ;;; look at the data, but only if there are obs
+        if (get_n_obs(out_arr) .gt. 0) then
+          summarize_array( \
+            this_fn+": summary of "+ sdate + "->" + out_dv_name + ":", \
+            sigdigs, out_arr)
+
+;           ;; View output in VERDI. Note:
+;           ;; * VERDI is sllooowwww when off LAN!     
+;           ;; * '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(out_fp_ncl, out_dv_name, False)
+        else
+          print(this_fn + ": ERROR?: no obs in " + out_fp_ncl)
+        end if ; (num(.not. ismissing(out_dv)) .gt. 0) then
+        delete(out_fh)
+      end if ; (first_day_this_month .eq. n_day)
+
+; ----------------------------------------------------------------------
+; teardown/cleanup
+; ----------------------------------------------------------------------
+
+      if ((.not. (isfilepresent(out_fp_cmaq))) .and. \
+          (.not. (isfilepresent(out_fp_ncl)))) then
+        print(this_fn+": ERROR: cannot find either monthly_CMAQ or monthly_NCL output")
+        return ; TODO: abend  
+      else if ((.not. (isfilepresent(out_fp_cmaq))) .and. \
+               (isfilepresent(out_fp_ncl))) then
+        ;; make file have the suffix we need :-( TODO: less system-dependently?
+        cmd = "mv "+out_fp_ncl+" "+out_fp_cmaq
+  ;       print(this_fn+": about to do: '"+cmd+"'") ; debugging
+        system(cmd)
+      else if ((isfilepresent(out_fp_cmaq)) .and. \
+               (isfilepresent(out_fp_ncl))) then
+        ;; remove file we don't need :-( TODO: less system-dependently?
+        cmd = "rm "+out_fp_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(out_fp_cmaq)
+      ; ASSERT: ((isfilepresent(out_fp_cmaq)) .and.
+      ;          (.not. (isfilepresent(out_fp_ncl))))
+
+    end do ; i_day = first_day_this_month, last_day_of_month
+  end do ; n_month = 1, 12
 
 end ; make_hourlies.ncl
Tip: Filter by directory path e.g. /media app.js to search for public/media/app.js.
Tip: Use camelCasing e.g. ProjME to search for ProjectModifiedEvent.java.
Tip: Filter by extension type e.g. /repo .js to search for all .js files in the /repo directory.
Tip: Separate your search with spaces e.g. /ssh pom.xml to search for src/ssh/pom.xml.
Tip: Use ↑ and ↓ arrow keys to navigate and return to view the file.
Tip: You can also navigate files with Ctrl+j (next) and Ctrl+k (previous) and view the file with Ctrl+o.
Tip: You can also navigate files with Alt+j (next) and Alt+k (previous) and view the file with Alt+o.