Source

CLM_CN_global_to_AQMEII-NA / reunit.ncl

;!/usr/bin/env ncl ; requires version >= 5.1.1
; Convert raw CLM-CN data from flux-rate to mass-rate prior to regridding,
; and exclude bogus "months":

; Per Eri Saikawa, input repeats:
; * there are only 12 "real" months
; * month 1==month 2
; * "month 13"=="month 14"
; Since NCL indexing is {C-style,0-based}, ignore bogus month indices 0,13

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

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

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

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

begin ; skip if copy/paste-ing to console

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

  this_fn = getenv("CALL_REUNIT_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"))
  CMAQ_radius = stringtofloat(getenv("CMAQ_RADIUS")) ; in m

  ;;; input: raw CLM-CN in flux-rate units (mgN/m^2)
  in_fp = getenv("CLM_CN_RAW_FP")
  in_datavar_name = getenv("CLM_CN_RAW_DATAVAR_NAME")
  in_datavar_longname = getenv("CLM_CN_RAW_DATAVAR_LONGNAME")
  in_datavar_unit = getenv("CLM_CN_RAW_DATAVAR_UNIT")
  ; in_datavar_na = stringtoint(getenv("CLM_CN_RAW_DATAVAR_NA"))
  in_x_dim_name = getenv("CLM_CN_RAW_X_COORDVAR_NAME")
  in_y_dim_name = getenv("CLM_CN_RAW_Y_COORDVAR_NAME")
  in_time_dim_name = getenv("CLM_CN_RAW_TIME_COORDVAR_NAME")

  ;;; output: same CLM-CN grid, but in mass-rate units (molN2O/s)
  out_fp = getenv("CLM_CN_REUNIT_FP")
  out_history_attr = getenv("CLM_CN_TARGET_FILE_ATTR_HISTORY")
  ;; datavar
  out_datavar_name = getenv("CLM_CN_REUNIT_DATAVAR_NAME")
  ;; IOAPI pads varattr=long_name to length=16 with trailing spaces
  out_datavar_long_name = getenv("CLM_CN_REUNIT_DATAVAR_LONG_NAME")
  ;; IOAPI pads varattr=units to length=16 with trailing spaces
  out_datavar_units = getenv("CLM_CN_REUNIT_DATAVAR_UNITS")
  ;; IOAPI pads varattr=var_desc to length=80 with trailing spaces
  out_datavar_var_desc = getenv("CLM_CN_REUNIT_DATAVAR_VAR_DESC")
  out_x_dim_name = getenv("CLM_CN_REUNIT_X_COORDVAR_NAME")
  out_y_dim_name = getenv("CLM_CN_REUNIT_Y_COORDVAR_NAME")
  out_time_dim_name = getenv("CLM_CN_REUNIT_TIME_COORDVAR_NAME")

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

;----------------------------------------------------------------------
; get input
;----------------------------------------------------------------------

  in_fh = addfile(in_fp, "r")
  in_dv = in_fh->$in_datavar_name$
  ;;; see note on input above
;  in_arr = (/ in_dv(1:12,:,:) /) ; minus metadata
  in_arr = in_dv(1:12,:,:)
  in_dims = dimsizes(in_arr)

; ; start debug---------------------------------------------------------
; printVarSummary(in_dv)
; ; Variable: in_dv
; ; Type: float
; ; Total Size: 774144 bytes
; ;             193536 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes:   [time | 14] x [lat | 96] x [lon | 144]
; ; Coordinates: 
; ;             time: [732892..733289]
; ;             lat: [-89.0625..89.0625]
; ;             lon: [1.25..358.75]
; ; Number Of Attributes: 2
; ;   long_name :   N2O emissions
; ;   units :       mgN/m2/month

; print(this_fn+": in_dims==")
; print(in_dims)
; Dimensions and sizes: [3]
; Coordinates: 
; (0)   12
; (1)   96
; (2)   144

; print(this_fn+": summarize_monthly_3d(in_arr)==")
; summarize_monthly_3d(in_arr)
; ; (0) month= 1: min=0.000e+00, avg=5.081e-01, max=4.511e+02
; ; (0) month= 2: min=0.000e+00, avg=5.596e-01, max=5.261e+02
; ; (0) month= 3: min=0.000e+00, avg=6.176e-01, max=4.793e+02
; ; (0) month= 4: min=0.000e+00, avg=6.738e-01, max=5.232e+02
; ; (0) month= 5: min=0.000e+00, avg=7.490e-01, max=6.628e+02
; ; (0) month= 6: min=0.000e+00, avg=1.098e+00, max=1.184e+03
; ; (0) month= 7: min=0.000e+00, avg=1.317e+00, max=1.012e+03
; ; (0) month= 8: min=0.000e+00, avg=1.519e+00, max=1.961e+03
; ; (0) month= 9: min=0.000e+00, avg=1.181e+00, max=5.836e+02
; ; (0) month=10: min=0.000e+00, avg=8.381e-01, max=5.904e+02
; ; (0) month=11: min=0.000e+00, avg=6.329e-01, max=5.078e+02
; ; (0) month=12: min=0.000e+00, avg=5.944e-01, max=9.043e+02

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

  ;----------------------------------------------------------------------
  ; get input areas
  ;----------------------------------------------------------------------

  ;;; get array of grid centers
;  lat_lon_arr = (/ in_arr(0,:,:) /) ; since only a template, can use any month
  ; but need to have named dimensions
  lat_lon_arr = in_dv(0,:,:)
  ;;; calculate gridcell areas in m^2
  areas_arr = lats_lons_sphere_area(lat_lon_arr, CMAQ_radius)

; ; start debug---------------------------------------------------------
; printVarSummary(areas_arr)
; ; Variable: areas_arr
; ; Type: float
; ; Total Size: 55296 bytes
; ;             13824 values
; ; Number of Dimensions: 2
; ; Dimensions and sizes:       [96] x [144]
; ; Coordinates: 
; ; Number Of Attributes: 1
; ;   _FillValue :      9.96921e+36
; ;   end debug---------------------------------------------------------

  ;----------------------------------------------------------------------
  ; calculate output data in mass-rate
  ;----------------------------------------------------------------------

  ;;; see
  ;;; Simulation-of-N2O-Production-and-Transport-in-the-US-Cornbelt-Compared-to-Tower-Measurements#wiki-input-processing-CLM-CNv3.5
  ;;; section="convert concentrations")

  ;;; molN2O   mgN      m^2        mo   gN    gN2O   molN2O
  ;;; ------ = ------ * -------- * -- * --- * ---- * ------
  ;;;      s   m^2 mo   gridcell   s    mgN   gN     gN2O

  ;;; Since the conversion is month-dependent, copy by month

  out_arr = new(dimsizes(in_arr), float)
  do i_month=0 , 11 ; remember: NCL indexing is {C-style,0-based}
    n_month = i_month + 1 ; 1 .. 12
    months_per_sec = 1.0 / tofloat(seconds_in_month(model_year, n_month))
    out_arr(i_month,:,:) = \
      tofloat(in_arr(i_month,:,:) * \        ; mgN/m^2/mo
      areas_arr * \                          ; m^2/gridcell
      1e-3 * \                               ; gN/mgN
      months_per_sec * \                     ; mo/s
      (molar_mass_n2o) / nitrogen_mass_n2o)  ; molN2O/gN
  end do ; iterate months

  ;;; Set the dimension/coordvar names/values
  out_arr!0 = out_time_dim_name
  out_arr&$out_time_dim_name$ = ispan(1,12,1)
  out_arr!1 = out_y_dim_name
  out_arr&$out_y_dim_name$ = in_dv&$in_y_dim_name$
  out_arr!2 = out_x_dim_name
  out_arr&$out_x_dim_name$ = in_dv&$in_x_dim_name$

; start debug---------------------------------------------------------
; printVarSummary(out_arr)
; Variable: out_arr
; Type: float
; Total Size: 663552 bytes
;             165888 values
; Number of Dimensions: 3
; Dimensions and sizes: [month | 12] x [lat | 96] x [lon | 144]
; Coordinates: 
;             month: [1..12]
;             lat: [-89.0625..89.0625]
;             lon: [1.25..358.75]
; Number Of Attributes: 1
;   _FillValue :        9.96921e+36

; print(this_fn+": summarize_monthly_3d(out_arr)==") ; in mol/s
; summarize_monthly_3d(out_arr)
; ; (0)   month= 1: min=0.000e+00, avg=1.627e+01, max=1.532e+04
; ; (0)   month= 2: min=0.000e+00, avg=1.921e+01, max=1.909e+04
; ; (0)   month= 3: min=0.000e+00, avg=1.980e+01, max=1.627e+04
; ; (0)   month= 4: min=0.000e+00, avg=2.167e+01, max=1.835e+04
; ; (0)   month= 5: min=0.000e+00, avg=2.098e+01, max=1.282e+04
; ; (0)   month= 6: min=0.000e+00, avg=3.048e+01, max=2.366e+04
; ; (0)   month= 7: min=0.000e+00, avg=3.655e+01, max=3.108e+04
; ; (0)   month= 8: min=0.000e+00, avg=4.269e+01, max=6.026e+04
; ; (0)   month= 9: min=0.000e+00, avg=3.543e+01, max=2.047e+04
; ; (0)   month=10: min=0.000e+00, avg=2.556e+01, max=2.004e+04
; ; (0)   month=11: min=0.000e+00, avg=2.068e+01, max=1.781e+04
; ; (0)   month=12: min=0.000e+00, avg=1.905e+01, max=3.070e+04
;   end debug---------------------------------------------------------
  
  ; ----------------------------------------------------------------------
  ; create output data, schema
  ; ----------------------------------------------------------------------

  ; 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}
  out_fh_att = True ; assign file/global attributes
  out_fh_att@HISTORY = out_history_attr
  out_fh_att@YEAR = model_year

  ;;; datavar and attributes
  ; 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(in_dims, float)
  ; ... the rest of the old datavar's metadata? Nope: temp_var must
  ; * have coordvars
  ; * have *the same* coordvars
;  copy_VarMeta(in_dv, temp_var)
  ; Since this is a non-file variable (at this point), we can delete attributes ...
  delete_VarAtts(temp_var, (/ "units", "long_name", "var_desc" /) )
  temp_var@units = out_datavar_units
  temp_var@long_name = out_datavar_long_name
  temp_var@var_desc = out_datavar_var_desc
  ; and copy the actual data :-)
;  temp_var = (/ out_arr /)
  ; ... and copy the coordvars
  temp_var = out_arr

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

; print(str_get_nl() + this_fn + ": printVarSummary(temp_var)==")
; printVarSummary(temp_var)
; ; Variable: temp_var
; ; Type: float
; ; Total Size: 663552 bytes
; ;             165888 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes:       [month | 12] x [lat | 96] x [lon | 144]
; ; Coordinates: 
; ;             month: [1..12]
; ;             lat: [-89.0625..89.0625]
; ;             lon: [1.25..358.75]
; ; Number Of Attributes: 4
; ;   var_desc :        Model species N2O                                                               
; ;   long_name :       N2O             
; ;   units :   moles/s         
; ;   _FillValue :      9.96921e+36

; print(str_get_nl() + this_fn + ": summarize_monthly_3d(temp_var)==")
; summarize_monthly_3d(temp_var)
; ; (0) month= 1: min=0.000e+00, avg=1.627e+01, max=1.532e+04
; ; (0) month= 2: min=0.000e+00, avg=1.921e+01, max=1.909e+04
; ; (0) month= 3: min=0.000e+00, avg=1.980e+01, max=1.627e+04
; ; (0) month= 4: min=0.000e+00, avg=2.167e+01, max=1.835e+04
; ; (0) month= 5: min=0.000e+00, avg=2.098e+01, max=1.282e+04
; ; (0) month= 6: min=0.000e+00, avg=3.048e+01, max=2.366e+04
; ; (0) month= 7: min=0.000e+00, avg=3.655e+01, max=3.108e+04
; ; (0) month= 8: min=0.000e+00, avg=4.269e+01, max=6.026e+04
; ; (0) month= 9: min=0.000e+00, avg=3.543e+01, max=2.047e+04
; ; (0) month=10: min=0.000e+00, avg=2.556e+01, max=2.004e+04
; ; (0) month=11: min=0.000e+00, avg=2.068e+01, max=1.781e+04
; ; (0) month=12: min=0.000e+00, avg=1.905e+01, max=3.070e+04

; ;   end debug---------------------------------------------------------
  
  ; ----------------------------------------------------------------------
  ; write output
  ; ----------------------------------------------------------------------

  out_fh_att@creation_date = systemfunc("date") ; TODO: non-system-dependent date
  out_fh = addfile(out_fp, "c")

  ;;; write global attributes
  fileattdef(out_fh, out_fh_att) ; copy file attributes   

  ;;; write datavar
  out_fh->$out_datavar_name$ = temp_var

; start debug---------------------------------------------------------
print(str_get_nl() + this_fn + ": printVarSummary(out_fh)==")
printVarSummary(out_fh)
; Variable: out_fh
; Type: file
; File path:                   .../2008PTONCLMCNN2O_reunit.nc
; Number of global attributes: 3
; Number of dimensions:        3
; Number of variables:         4

print(str_get_nl() + this_fn + ": printFileVarSummary(out_fh, out_datavar_name)==")
printFileVarSummary(out_fh, out_datavar_name)
; Variable: N2O
; Type: float
; Total Size: 663552 bytes
;             165888 values
; Number of Dimensions: 3
; Dimensions and sizes:	[month | 12] x [lat | 96] x [lon | 144]
; Coordinates: 
;             month: [1..12]
;             lat: [-89.0625..89.0625]
;             lon: [1.25..358.75]
; Number of Attributes: 4
;   var_desc :	Model species N2O                                                               
;   long_name :	N2O             
;   units :	moles/s         
;   _FillValue :	9.96921e+36

print(str_get_nl() + this_fn + ": summarize_monthly_3d(out_fh->out_datavar_name)==")
summarize_monthly_3d(out_fh->$out_datavar_name$)
; (0)   month= 1: min=0.000e+00, avg=1.627e+01, max=1.532e+04
; (0)   month= 2: min=0.000e+00, avg=1.921e+01, max=1.909e+04
; (0)   month= 3: min=0.000e+00, avg=1.980e+01, max=1.627e+04
; (0)   month= 4: min=0.000e+00, avg=2.167e+01, max=1.835e+04
; (0)   month= 5: min=0.000e+00, avg=2.098e+01, max=1.282e+04
; (0)   month= 6: min=0.000e+00, avg=3.048e+01, max=2.366e+04
; (0)   month= 7: min=0.000e+00, avg=3.655e+01, max=3.108e+04
; (0)   month= 8: min=0.000e+00, avg=4.269e+01, max=6.026e+04
; (0)   month= 9: min=0.000e+00, avg=3.543e+01, max=2.047e+04
; (0)   month=10: min=0.000e+00, avg=2.556e+01, max=2.004e+04
; (0)   month=11: min=0.000e+00, avg=2.068e+01, max=1.781e+04
; (0)   month=12: min=0.000e+00, avg=1.905e+01, max=3.070e+04
;   end debug---------------------------------------------------------

end ; reunit.ncl