Source

GEIA_regrid / retemp_reunit.ncl

Full commit
;!/usr/bin/env ncl ; requires version >= 5.1.1
;----------------------------------------------------------------------
; Convert regridded GEIA data to the units and temporality required for CMAQ.

;----------------------------------------------------------------------

;----------------------------------------------------------------------
; 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 int2flt,  *_VarAtts

;----------------------------------------------------------------------
; functions
;----------------------------------------------------------------------

; can't call `getenv` before `begin`?
load "./summarize.ncl"
load "./time.ncl"

;----------------------------------------------------------------------
; code
;----------------------------------------------------------------------

begin ; skip if copy/paste-ing to console

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

  this_fn = getenv("CALL_RETEMP_FN") ; for debugging

  ;; generic model constants
  model_year = stringtoint(getenv("MODEL_YEAR"))
  molar_mass_N2O = stringtofloat(getenv("MOLAR_MASS_N2O"))
  nitrogen_mass_N2O = stringtofloat(getenv("NITROGEN_MASS_N2O"))
  grams_per_ton = stringtofloat(getenv("GRAMS_PER_TON"))

  ;; regridded output from R script is our input
  regrid_var_name = getenv("GEIA_REGRID_DATAVAR_NAME")
  regrid_fp = getenv("GEIA_REGRID_FP")    ; path

  ;; create data container files from extents/template file
  template_fp = getenv("TEMPLATE_INPUT_FP")
  template_fp_var_name = getenv("TEMPLATE_VAR_NAME")
  ; driver must
  ; * handle fact that NCL {does not support, won't create} *.ncf :-(
  ; * ensure output_fp_template* do not exist!
  output_fp_template_CMAQ = getenv("GEIA_TARGET_FP_TEMPLATE_CMAQ")
  output_fp_template_NCL = getenv("GEIA_TARGET_FP_TEMPLATE_NCL")
  ; replace output_fp_template_str with <year/><month/> to create "real" filepaths
  output_fp_template_str = getenv("GEIA_TARGET_FN_TEMPLATE_STR")
  output_history_attr = getenv("GEIA_TARGET_FILE_ATTR_HISTORY")
  ; datavar metadata
  output_fp_template_var_name = getenv("GEIA_TARGET_VAR_NAME")
  output_fp_template_var_attr_long_name = getenv("GEIA_TARGET_VAR_ATTR_LONG_NAME")
  output_fp_template_var_attr_units = getenv("GEIA_TARGET_VAR_ATTR_UNITS")
  output_fp_template_var_attr_var_desc = getenv("GEIA_TARGET_VAR_ATTR_VAR_DESC")

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

  ;;; setup template file for yearly output
  template_fh = addfile(template_fp, "r")                  ; file handle
  output_fp_template_fh = addfile(output_fp_template_NCL, "c")
;  print(output_fp_template_NCL) ; debugging

  ;; 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.
  output_fp_template_fh_att = True ; assign file/global attributes
  output_fp_template_fh_att@HISTORY = output_history_attr
  output_fp_template_fh_att@YEAR = model_year

  ;; 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 old datavar), and ...
  temp_var = new( \
    dimsizes(template_fh->$template_fp_var_name$), \
    float, "No_FillValue") ; since no _FillValue in original
  ; ... the rest of the old datavar's metadata.
  copy_VarMeta(template_fh->$template_fp_var_name$, temp_var)
  ; Since this is a non-file variable (at this point), we can delete attributes ...
  delete_VarAtts(temp_var, (/ "long_name", "var_desc" /) )
  temp_var@long_name = output_fp_template_var_attr_long_name
  temp_var@var_desc = output_fp_template_var_attr_var_desc
  ; ... and then write the "new datavar" to the file
  output_fp_template_fh->$output_fp_template_var_name$ = temp_var
  delete(temp_var)

; ; start debug---------------------------------------------------------

; printVarSummary(output_fp_template_fh)
; ; Variable: output_fp_template_fh
; ; Type: file
; ; File path:    /home/rtd/code/regridding/GEIA_regrid/emis_mole_N2O_######_12US1_cmaq_cb05_soa_2008ab_08c.ncf.nc
; ; Number of global attributes:   0
; ; Number of dimensions:  4
; ; Number of variables:   1

; printFileVarSummary(output_fp_template_fh, output_fp_template_var_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: not a coordinate variable
; ;             LAY: not a coordinate variable
; ;             ROW: not a coordinate variable
; ;             COL: not a coordinate variable
; ; Number of Attributes: 3
; ;   var_desc :  Model species N2O                                                               
; ;   long_name : N2O             
; ;   units :     moles/s         

; ;   end debug---------------------------------------------------------

  ;;; reunit/retemp the previously-regridded output
  regrid_fh = addfile(regrid_fp, "r") ; handle
  regrid_var = regrid_fh->$regrid_var_name$

  ;; calculate conversion coefficient
  ; input GEIA gridcell/timestep N2O (in (metric) tonN/yr/gridcell)
  ;   / sec/yr  ; varies with year
  ;   * gN2O/gN
  ;   * gN/tonN ; aka g/ton
  ;   = output gridcell/timestep N2O (in moles/s), where
  ; input timestep==yearly
  ; output timestep==hourly
  reunit_arr = regrid_var(:,:) * molar_mass_N2O * tofloat(grams_per_ton) /\
    (tofloat(seconds_in_year(model_year)) * nitrogen_mass_N2O)

; ; start debug---------------------------------------------------------

; printVarSummary(regrid_var)
; ; Variable: regrid_var
; ; Type: float
; ; Total Size: 548964 bytes
; ;             137241 values
; ; Number of Dimensions: 2
; ; Dimensions and sizes: [ROW | 299] x [COL | 459]
; ; Coordinates: 
; ;             ROW: [1854000..-1722000]
; ;             COL: [-2550000..2946000]
; ; Number Of Attributes: 8
; ;   _FillValue :        -999
; ;   units :     ton N2O-N/yr
; ;   missing_value :     -999
; ;   long_name : N2O emissions
; ;   projection :        +proj=lcc +lat_1=33 +lat_2=45 +lat_0=40 +lon_0=-97 +a=6370000 +b=6370000
; ;   projection_format : PROJ.4
; ;   min :       0.83728
; ;   max :       522.6937746382763

; printVarSummary(reunit_arr)
; ; Variable: reunit_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 :	-999

; ;   end debug---------------------------------------------------------

  ;;; write output to annual file

  ;; copy template file to annual file
  output_fp_NCL = \
    str_sub_str(output_fp_template_NCL, output_fp_template_str, \
      sprinti("%4i", model_year))
;      sprinti("%4i", model_year)+sprinti("%0.2i", month))
  output_fp_CMAQ = \
    str_sub_str(output_fp_template_CMAQ, output_fp_template_str, \
      sprinti("%4i", model_year))
  output_fp_NCL_fh = addfile(output_fp_NCL, "c") ; handle

  ;; write global attributes to annual file
  output_fp_fh_att = output_fp_template_fh_att
  ; TODO: non-system-dependent date
  output_fp_fh_att@creation_date = systemfunc("date")
  fileattdef( \
    output_fp_NCL_fh, \
    output_fp_fh_att) ; copy file attributes   

  ;; write datavar definition to annual file
;    copy_VarMeta( \
;      output_fp_template_fh->$output_fp_template_var_name$, \
;      output_fp_NCL_fh->$output_fp_template_var_name$)
  ; gotta use the temporary datavar trick to copy (again)
  temp_var = new( \
    dimsizes(output_fp_template_fh->$output_fp_template_var_name$), \
    float, "No_FillValue") ; since no _FillValue in original
  copy_VarMeta(output_fp_template_fh->$output_fp_template_var_name$, temp_var)
  ; ... and then write the "new datavar" to the file
  output_fp_NCL_fh->$output_fp_template_var_name$ = temp_var
;  printVarSummary(output_fp_NCL_fh) ; debugging

  ;; write 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 | :)
  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_fp_NCL_fh->$output_fp_template_var_name$(i_hour,0,:,:) = (/ reunit_arr(:,:) /)
  end do
  
  ; rename the {reunit,retemp}ed annual file. TODO: make non-system-dependent.
  cmd = "mv "+output_fp_NCL+" "+output_fp_CMAQ
;  print(this_fn+": about to do: '"+cmd+"'") ; debugging
  system(cmd)

; ; start debug---------------------------------------------------------

;   output_fp_NCL_fh = addfile(output_fp_NCL, "r") ; handle
;   printVarSummary(output_fp_NCL_fh)
;   printVarSummary(output_fp_NCL_fh->$output_fp_template_var_name$)
;   print(str_get_nl() + this_fn +": about to do: '" +\
;     "summarize_annual_4d(output_fp_NCL_fh->" +\
;     output_fp_template_var_name + "'")
;   summarize_annual_4d(output_fp_NCL_fh->$output_fp_template_var_name$)

; ;   end debug---------------------------------------------------------

  ;;; cleanup/teardown

end ; retemp_reunit.ncl