Source

CLM_CN_global_to_AQMEII-NA / reunit.ncl

;!/usr/bin/env ncl ; requires version >= 5.1.1
;----------------------------------------------------------------------
; Convert regridded CLM-CN 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`?
; can't even load from a string variable?
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"))
  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_var_name = getenv("CLM_CN_REGRID_DATAVAR_NAME")
  regrid_fp = getenv("CLM_CN_REGRID_FP")    ; path

  ;; load matrix of gridcell map scale factors (MSFs) from GRIDCRO2D_<date/>
  msf_fp = getenv("MSF_FP")    ; path to container netCDF file
  msf_var_name = getenv("MSF_VAR_NAME")

  ;; 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("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")
  output_history_attr = getenv("CLM_CN_TARGET_FILE_ATTR_HISTORY")
  ; datavar metadata
  output_fp_template_var_name = getenv("CLM_CN_TARGET_VAR_NAME")
  output_fp_template_var_attr_long_name = getenv("CLM_CN_TARGET_VAR_ATTR_LONG_NAME")
  output_fp_template_var_attr_units = getenv("CLM_CN_TARGET_VAR_ATTR_UNITS")
  output_fp_template_var_attr_var_desc = getenv("CLM_CN_TARGET_VAR_ATTR_VAR_DESC")

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

  ;;; setup template file for monthly output
  template_fh = addfile(template_fp, "r")                  ; file handle
  output_fp_template_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)
  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
; do once for monthlies
;  fileattdef(output_fp_template_fh, output_fp_template_fh_att) ; copy file attributes   

  ;; 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:  ./emis_mole_N2O_######_12US1_cmaq_cb05_soa_2008ab_08c.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 (file variable)
; ; 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: 
; ; Number Of Attributes: 3
; ;   var_desc :        Model species N2O                                                               
; ;   long_name :       N2O             
; ;   units :   moles/s         

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

  ;;; delete bogus "months" 1,14 in regrid output

  regrid_fh = addfile(regrid_fp,"r") ; handle
  regrid_var = regrid_fh->$regrid_var_name$
  ; keep non-bogus months: remember NCL indexing is {C-style,0-based}
  regrid_arr = regrid_var(1:12,:,:)

; ; start debug---------------------------------------------------------
; printVarSummary(regrid_arr) 
; ; Variable: regrid_arr
; ; Type: float
; ; Total Size: 6587568 bytes
; ;             1646892 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes: [time | 12] x [ROW | 299] x [COL | 459]
; ; Coordinates: 
; ;             time: [2..13]
; ;             ROW: [1854000..-1722000]
; ;             COL: [-2550000..2946000]
; ; Number Of Attributes: 8
; ;   max :       <ARRAY of 14 elements>
; ;   min :       <ARRAY of 14 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 :     mgN/m2/month
; ;   _FillValue :        -3.4e+38

; summarize_monthly_3d(regrid_arr) 
; ; (0)   month= 1: min=0.000e+00, avg=3.540e-01, max=5.116e+00
; ; (0)   month= 2: min=0.000e+00, avg=3.689e-01, max=7.540e+00
; ; (0)   month= 3: min=0.000e+00, avg=5.577e-01, max=3.397e+01
; ; (0)   month= 4: min=0.000e+00, avg=1.086e+00, max=3.558e+01
; ; (0)   month= 5: min=0.000e+00, avg=1.728e+00, max=2.909e+01
; ; (0)   month= 6: min=0.000e+00, avg=3.068e+00, max=1.500e+02
; ; (0)   month= 7: min=0.000e+00, avg=3.193e+00, max=1.742e+02
; ; (0)   month= 8: min=0.000e+00, avg=2.962e+00, max=4.402e+01
; ; (0)   month= 9: min=0.000e+00, avg=2.745e+00, max=5.323e+01
; ; (0)   month=10: min=0.000e+00, avg=1.724e+00, max=2.571e+01
; ; (0)   month=11: min=0.000e+00, avg=8.248e-01, max=1.455e+01
; ; (0)   month=12: min=0.000e+00, avg=5.972e-01, max=1.520e+01

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

  ;;; reunit and retemporalize cleaned 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")

  ;; load gridcell areas = nominal area * map scale factor (MSF)
  msf_fh = addfile(msf_fp,"r") ; handle
  msf_var = msf_fh->$msf_var_name$
  ; ignore the unitary layer and timestep, multiply by nominal area
  area_arr = nominal_gridcell_area * msf_var(0,0,:,:)

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

; ; printFileVarSummary(msf_fh, msf_var_name)
; printVarSummary(msf_var)
; ; ; Variable: msf_var
; ; ; Type: float
; ; ; Total Size: 548964 bytes
; ; ;             137241 values
; ; ; Number of Dimensions: 4
; ; ; Dimensions and sizes: [TSTEP | 1] x [LAY | 1] x [ROW | 299] x [COL | 459]
; ; ; Coordinates: 
; ; ; Number Of Attributes: 3
; ; ;   long_name : MSFX2           
; ; ;   units :     (M/M)**2        
; ; ;   var_desc :  squared map-scale factor (CROSS)

; printVarSummary(area_arr) 
; ; Variable: area_arr
; ; Type: float
; ; Total Size: 548964 bytes
; ;             137241 values
; ; Number of Dimensions: 2
; ; Dimensions and sizes:       [299] x [459]
; ; Coordinates: 

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

  ;; process gridcells by month

  ; copy structure of regrid_arr, create new data container (anality?)
  reunit_arr = new(dimsizes(regrid_arr), float)
  do i_month=0 , 11 ; remember: NCL indexing is {C-style,0-based}
    reunit_arr(i_month,:,:) = \
      tofloat(regrid_arr(i_month,:,:) * area_arr * molar_mass_n2o) /\
      (tofloat(seconds_in_month(2008,i_month+1)) * nitrogen_mass_n2o)
  end do ; iterate months

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

; summarize_monthly_3d(reunit_arr) ; moles/s/gridcell over all gridcells
; ; month= 1: min=0.000e+00, avg=3.028e+01, max=4.477e+02
; ; month= 2: min=0.000e+00, avg=3.377e+01, max=6.986e+02
; ; month= 3: min=0.000e+00, avg=4.766e+01, max=2.944e+03
; ; month= 4: min=0.000e+00, avg=9.568e+01, max=3.186e+03
; ; month= 5: min=0.000e+00, avg=1.473e+02, max=2.515e+03
; ; month= 6: min=0.000e+00, avg=2.695e+02, max=1.328e+04
; ; month= 7: min=0.000e+00, avg=2.720e+02, max=1.501e+04
; ; month= 8: min=0.000e+00, avg=2.523e+02, max=3.806e+03
; ; month= 9: min=0.000e+00, avg=2.414e+02, max=4.597e+03
; ; month=10: min=0.000e+00, avg=1.469e+02, max=2.157e+03
; ; month=11: min=0.000e+00, avg=7.269e+01, max=1.303e+03
; ; month=12: min=0.000e+00, avg=5.101e+01, max=1.317e+03

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

  ;;; write output to monthly files
  do i_month=0 , 11 ; remember: NCL indexing is {C-style,0-based}
    ;; copy template file to monthly file
    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", month))
    output_fp_monthly_CMAQ = \
      str_sub_str(output_fp_template_CMAQ, output_fp_template_str, \
        sprinti("%4i", model_year)+sprinti("%0.2i", month))
    output_fp_monthly_NCL_fh = \
      addfile(output_fp_monthly_NCL, "c") ; handle

    ;; write global attributes to monthly file
    output_fp_monthly_fh_att = output_fp_template_fh_att
    ; TODO: non-system-dependent date
    output_fp_monthly_fh_att@creation_date = systemfunc("date")
    output_fp_monthly_fh_att@MONTH = month
    fileattdef( \
      output_fp_monthly_NCL_fh, \
      output_fp_monthly_fh_att) ; copy file attributes   

    ;; write datavar definition to monthly file
;     copy_VarMeta( \
;       output_fp_template_fh->$output_fp_template_var_name$, \
;       output_fp_monthly_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_monthly_NCL_fh->$output_fp_template_var_name$ = temp_var

    ;; 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
      output_fp_monthly_NCL_fh->$output_fp_template_var_name$(i_hour,0,:,:) = (/ reunit_arr(i_month,:,:) /)
    end do

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

; print(str_get_nl() + str_get_nl() + \
;   this_fn+": output from month="+sprinti("%0.2i", month)+":")
; ; printVarSummary(output_fp_monthly_NCL_fh)
; printFileVarSummary(output_fp_monthly_NCL_fh, output_fp_template_var_name)
; summarize_monthly_4d(\
;   output_fp_monthly_NCL_fh->$output_fp_template_var_name$, month)

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

    ; 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_monthly_to_hourly.ncl