Source

CLM_CN_global_to_AQMEII-NA / retemp.ncl

Full commit
;!/usr/bin/env ncl ; requires version >= 5.1.1
; Convert regridded CLM-CN data to the 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 

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

; can't call `getenv` before `begin`?
; can't even load from a string variable?
load "$SUMMARIZE_FUNCS_FP"
load "$TIME_FUNCS_FP"

;----------------------------------------------------------------------
; 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"))
  nominal_gridcell_area = stringtofloat(getenv("NOMINAL_GRIDCELL_AREA"))
  molar_mass_n2o = stringtofloat(getenv("MOLAR_MASS_N2O"))
  nitrogen_mass_n2o = stringtofloat(getenv("NITROGEN_MASS_N2O"))
  seconds_per_hour = stringtofloat(getenv("SECONDS_PER_HOUR"))

  ;; regridded output from R script is our input
  regrid_fp = getenv("CLM_CN_REGRID_FP")    ; path
  regrid_datavar_name = getenv("CLM_CN_REGRID_DATAVAR_NAME")
  regrid_coordvar_x_name = getenv("CLM_CN_REGRID_X_COORDVAR_NAME")
  regrid_coordvar_y_name = getenv("CLM_CN_REGRID_Y_COORDVAR_NAME")
  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")
  ; note template has no coordvars

  ; 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("CLM_CN_TARGET_FP_TEMPLATE_CMAQ")
  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")
  ; 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")

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

; system("cp ~/code/regridding/CLM_CN_global_to_AQMEII-NA/2008PTONCLMCNN2O_reunit.nc ./")                                   
; system("cp ~/code/regridding/CLM_CN_global_to_AQMEII-NA/2008PTONCLMCNN2O_reunit_regrid.nc ./")

;----------------------------------------------------------------------
; create template output
;----------------------------------------------------------------------

  ;;; 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
  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")

  ;;; write output to monthly files
  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" + \
      str_get_nl() + str_get_tab() + 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
    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 | :)
    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,:,:) = (/ regrid_dv(i_month,:,:) /)
    end do

    ; rename the monthly file. TODO: make non-system-dependent.
    cmd = "mv "+output_fp_monthly_NCL+" "+output_fp_monthly_CMAQ
;    print(this_fn+": about to do: '"+cmd+"'") ; debugging
    system(cmd)
;    delete(\
;      (/ cmd, output_fp_monthly_CMAQ, output_fp_monthly_NCL, \
;        output_fp_monthly_NCL_fh /) )
; > 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, \
         output_fp_monthly_NCL_fh /] )

  end do ; iterate months
  
  ;;; cleanup/teardown

end ; retemp.ncl