Source

GFED-3.1_global_to_AQMEII-NA / check_conservation.ncl

Full commit
;!/usr/bin/env ncl ; requires version >= ???
;----------------------------------------------------------------------
; Check conservation of mass from GFED global monthlies to regridded (by me) monthlies to regridded hourlies, using input units.

;----------------------------------------------------------------------
; libraries
;----------------------------------------------------------------------

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"    ; all built-ins?
;load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl" ; for stat_dispersion

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 several

;----------------------------------------------------------------------
; constants
;----------------------------------------------------------------------

;;; model constants

this_fn = getenv("CHECK_CONSERVATION_FN") ; for debugging
model_year = stringtoint(getenv("MODEL_YEAR"))
hours_per_day =  stringtoint(getenv("HOURS_PER_DAY"))
seconds_per_hour = stringtoint(getenv("SECONDS_PER_HOUR"))

;;; global monthlies (as molarized netCDF, not the rawest GFED ASCII)
global_monthly_dv_name = getenv("GFED_GLOBAL_MONTHLY_CONV_DATAVAR_NAME")
global_monthly_fp = getenv("GFED_GLOBAL_MONTHLY_CONV_FP")

;;; regridded monthlies
regrid_monthly_dv_name = getenv("GFED_REGRID_MONTHLY_DATAVAR_NAME")
regrid_monthly_fp = getenv("GFED_REGRID_MONTHLY_FP")

;;; regridded hourlies
hourly_dv_name = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME")
hourly_template_string = getenv("GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING")
hourly_fp_template = getenv("GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL")
hourly_arr_dims_row_n = stringtoint(getenv("GFED_REGRID_N_ROW"))
hourly_arr_dims_col_n = stringtoint(getenv("GFED_REGRID_N_COL"))
hourly_arr_dims_hour_n = hours_in_year(model_year)
hourly_arr_dims = \ ; TODO: check whether order affects performance
  (/ hourly_arr_dims_hour_n, hourly_arr_dims_row_n, hourly_arr_dims_col_n /)
i_hours_in_year = 0 ; will count to hours_in_year(model_year)-1

; conservation report constants
sigdigs = stringtoint(getenv("OUTPUT_SIGNIFICANT_DIGITS"))
report_field_width = stringtoint(getenv("GFED_CONSERV_REPORT_FIELD_WIDTH"))
report_col_sep = getenv("GFED_CONSERV_REPORT_COLUMN_SEPARATOR")
report_title = getenv("GFED_CONSERV_REPORT_TITLE")
report_subtitle = getenv("GFED_CONSERV_REPORT_SUBTITLE")
report_flt_format = getenv("GFED_CONSERV_REPORT_FLOAT_FORMAT")
report_int_format = getenv("GFED_CONSERV_REPORT_INT_FORMAT")

;----------------------------------------------------------------------
; functions
;----------------------------------------------------------------------
  
;----------------------------------------------------------------------
; code
;----------------------------------------------------------------------

begin ; skip if copy/paste-ing to console

;----------------------------------------------------------------------
; payload
;----------------------------------------------------------------------

;----------------------------------------------------------------------
; process datasets
;----------------------------------------------------------------------

;----------------------------------------------------------------------
; process global monthlies
;----------------------------------------------------------------------

  global_monthly_fh = addfile(global_monthly_fp, "r")
  global_monthly_dv = global_monthly_fh->$global_monthly_dv_name$
  global_monthly_arr = (/ global_monthly_dv /)
  global_monthly_sum = sum(global_monthly_arr)
  global_monthly_stats = stat_dispersion(global_monthly_arr, False) ; False == don't stdout
  global_monthly_cells = global_monthly_stats(18) ; "Total number of input values"
  global_monthly_nNA = global_monthly_stats(20)   ; "Number of missing values"
  global_monthly_obs = global_monthly_cells - global_monthly_nNA

  print(\
    str_get_nl()                                            + \
    "GFED monthly N2O (mol) over globe from " + global_monthly_fp + ":" + str_get_nl()  + \
    "|cells|=" + global_monthly_cells                       + str_get_nl() + \
    "|obs|  =" + global_monthly_obs                         + str_get_nl() + \
    "min    =" + sprintf("%7.3e", global_monthly_stats(2))  + str_get_nl() + \
    "q1     =" + sprintf("%7.3e", global_monthly_stats(6))  + str_get_nl() + \
    "med    =" + sprintf("%7.3e", global_monthly_stats(8))  + str_get_nl() + \
    "mean   =" + sprintf("%7.3e", global_monthly_stats(0))  + str_get_nl() + \
    "q3     =" + sprintf("%7.3e", global_monthly_stats(10)) + str_get_nl() + \
    "max    =" + sprintf("%7.3e", global_monthly_stats(14)) + str_get_nl() + \
    "sum    =" + sprintf("%7.3e", global_monthly_sum)       + str_get_nl() + \
    str_get_nl() )

;----------------------------------------------------------------------
; process regridded monthlies
;----------------------------------------------------------------------

  regrid_monthly_fh = addfile(regrid_monthly_fp, "r")
  regrid_monthly_dv = regrid_monthly_fh->$regrid_monthly_dv_name$
  regrid_monthly_arr = (/ regrid_monthly_dv /)
  regrid_monthly_sum = sum(regrid_monthly_arr)
  regrid_monthly_stats = stat_dispersion(regrid_monthly_arr, False)
  regrid_monthly_cells = regrid_monthly_stats(18) ; "Total number of input values"
  regrid_monthly_nNA = regrid_monthly_stats(20)   ; "Number of missing values"
  regrid_monthly_obs = regrid_monthly_cells - regrid_monthly_nNA

  print(\
    str_get_nl()                                            + \
    "monthly N2O (mol) over AQMEII-NA from " + regrid_monthly_fp + ":" + str_get_nl() + \
    "|cells|=" + regrid_monthly_cells                       + str_get_nl() + \
    "|obs|  =" + regrid_monthly_obs                         + str_get_nl() + \
    "min    =" + sprintf("%7.3e", regrid_monthly_stats(2))  + str_get_nl() + \
    "q1     =" + sprintf("%7.3e", regrid_monthly_stats(6))  + str_get_nl() + \
    "med    =" + sprintf("%7.3e", regrid_monthly_stats(8))  + str_get_nl() + \
    "mean   =" + sprintf("%7.3e", regrid_monthly_stats(0))  + str_get_nl() + \
    "q3     =" + sprintf("%7.3e", regrid_monthly_stats(10)) + str_get_nl() + \
    "max    =" + sprintf("%7.3e", regrid_monthly_stats(14)) + str_get_nl() + \
    "sum    =" + sprintf("%7.3e", regrid_monthly_sum)       + str_get_nl() + \
    str_get_nl() )

  print(str_get_nl() + "AQMEII-NA monthly N2O/global monthly N2O==" + \
    sprintf("%7.3e", (regrid_monthly_sum / global_monthly_sum)))
  
;----------------------------------------------------------------------
; process (regridded) hourlies
;----------------------------------------------------------------------

  ;; 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
    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
      hourly_fp = get_daily_emis_fp(\
        n_month, n_day, hourly_fp_template, hourly_template_string)
      hourly_fh = addfile(hourly_fp, "r")
      hourly_dv = hourly_fh->$hourly_dv_name$
      do n_hour = 1, hours_per_day ; note the files go to hours_per_day+1
        ; TODO: show progress
        i_hour = n_hour - 1
        hourly_arr(i_hours_in_year,:,:) = (/ hourly_dv(i_hour,0,:,:) /)
        i_hours_in_year = i_hours_in_year + 1
      end do ; iterating hours

      ; just to avoid spurious errors to console
      delete(hourly_fh)
      delete(hourly_dv)

    end do ; iterating days
  end do ; iterating months

  ; start debugging---------------------------------------------------
  print("hourly_arr==")
  printVarSummary(hourly_arr)
  ;   end debugging---------------------------------------------------

  ;; convert mol/s -> mol
  hourly_arr = hourly_arr * seconds_per_hour
  hourly_sum = sum(hourly_arr)
  ;;; stat_dispersion(...) takes too long, so do that after the simple report
  hourly_nNA = num(ismissing(hourly_arr)) ; |missing values|

;----------------------------------------------------------------------
; generate report
;----------------------------------------------------------------------

  ;;; conservation report, with unit conversion
  print(str_get_nl()+report_title+"")
;  print(report_subtitle) ; prints variable information!
  print(report_subtitle+"")

  print(\  ; header1
    pad_blanks_to_length("input", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("input", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("output", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("output", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("", -report_field_width) )

  print(\  ; header2
    pad_blanks_to_length("global", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("NAs", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("AQMEII-NA", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("NAs", -report_field_width) + report_col_sep +\
    pad_blanks_to_length("out/in", -report_field_width) )

  print(\  ; content
    sprintf(report_flt_format, global_monthly_sum)        + report_col_sep +\
    sprinti(report_int_format, toint(global_monthly_nNA)) + report_col_sep +\
    sprintf(report_flt_format, hourly_sum)                + report_col_sep +\
    sprinti(report_int_format, toint(hourly_nNA))         + report_col_sep +\
    sprintf(report_flt_format, hourly_sum/global_monthly_sum) )

  print(str_get_nl())

  ; start progress----------------------------------------------------
  ; TODO: better progress control
  print(str_get_nl() + \
    this_fn+": start summarizing aggregate of hourlies @ " + \
    systemfunc("date") + str_get_nl() + \
    "will take awhile! (15-60 min on terrae interactive)" + str_get_nl() )
  ;   end progress----------------------------------------------------

  hourly_stats = stat_dispersion(hourly_arr, False)
;  hourly_nNA = hourly_stats(20) ; |missing values|

  ; start progress----------------------------------------------------
  print(str_get_nl() + \
    this_fn+": ended summarizing aggregate of hourlies @ " + \
    systemfunc("date") + str_get_nl() )
  ;   end progress----------------------------------------------------

  print(\
    str_get_nl()                                        + \
    "AQMEII-NA N2O (mol) from output hourly emissions:" + str_get_nl() + \
    "|cells|=" + hourly_stats(18)                       + str_get_nl() + \
    "|obs|  =" + hourly_stats(19)                       + str_get_nl() + \
    "min    =" + sprintf("%7.3e", hourly_stats(2))      + str_get_nl() + \
    "q1     =" + sprintf("%7.3e", hourly_stats(6))      + str_get_nl() + \
    "med    =" + sprintf("%7.3e", hourly_stats(8))      + str_get_nl() + \
    "mean   =" + sprintf("%7.3e", hourly_stats(0))      + str_get_nl() + \
    "q3     =" + sprintf("%7.3e", hourly_stats(10))     + str_get_nl() + \
    "max    =" + sprintf("%7.3e", hourly_stats(14))     + str_get_nl() + \
    "sum    =" + sprintf("%7.3e", hourly_sum)           + str_get_nl() + \
    str_get_nl() )

end ; check_conservation.ncl