1. Tom Roche
  2. GFED-3.1_global_to_AQMEII-NA

Source

GFED-3.1_global_to_AQMEII-NA / make_hourlies.ncl

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
;!/usr/bin/env ncl ; requires version >= 5.1.1 for str_sub_str
;----------------------------------------------------------------------
; Multiply GFED monthly emissions by its various fractions to create 3-hourly emissions,
; and from those write CMAQ-suitable hourly emissions files, following instructions @
; ftp://gfed3:dailyandhourly@zea.ess.uci.edu/GFEDv3.1/Readme.pdf

;----------------------------------------------------------------------
; 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 
load "$WORK_DIR/get_daily_emis_fp.ncl" ; for that function
;load "$WORK_DIR/string.ncl" ; for right_pad_blanks
load "$WORK_DIR/time.ncl" ; for first_day_of_month, last_day_of_month

;----------------------------------------------------------------------
; globals
; define here to make available in function definitions
;----------------------------------------------------------------------

this_fn = getenv("MAKE_HOURLIES_FN") ; for debugging

;;; generic model constants
model_year = stringtoint(getenv("MODEL_YEAR"))
hours_per_day = stringtoint(getenv("HOURS_PER_DAY"))
seconds_per_hour = stringtofloat(getenv("SECONDS_PER_HOUR"))
molar_mass_N2O = stringtofloat(getenv("molar_mass_N2O"))

;;; inputs

;; from GFED (before regridding)
monthly_emis_fp = getenv("GFED_REGRID_MONTHLY_FP")
monthly_emis_dv_name = getenv("GFED_REGRID_MONTHLY_DATAVAR_NAME")
daily_frac_fp = getenv("GFED_REGRID_DAILY_FRAC_FP")
daily_frac_dv_name = getenv("GFED_REGRID_DAILY_FRAC_DATAVAR_NAME")
hour3ly_frac_fp = getenv("GFED_REGRID_3HOURLY_FRAC_FP")
hour3ly_frac_dv_name = getenv("GFED_REGRID_3HOURLY_FRAC_DATAVAR_NAME")

;; AQMEII-NA areas: thought I needed, but am leaving them in anyway
out_areas_fp = getenv("GFED_AREAS_AQMEIINA_FP")
out_areas_dv_name = getenv("GFED_AREAS_AQMEIINA_DATAVAR_NAME")

;;; outputs

;; file
; templates
day_template_string = getenv("GFED_REGRID_HOURLY_EMIS_TEMPLATE_STRING")
out_emis_fp_template_cmaq = getenv("GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_CMAQ")
out_emis_fp_template_ncl = getenv("GFED_REGRID_HOURLY_EMIS_FP_TEMPLATE_NCL")
; attributes: TODO: port function get_output_file_attributes
out_attr_file_history = getenv("GFED_REGRID_HOURLY_EMIS_HISTORY")

;; datavar
out_dv_name = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_NAME")
; dimensions
out_tstep_dim_n = stringtoint(getenv("GFED_REGRID_HOURLY_EMIS_N_TSTEP"))
out_lay_dim_n = stringtoint(getenv("GFED_REGRID_HOURLY_EMIS_N_LAYER"))
; attributes
out_dv_attr_units = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_UNITS")
out_dv_attr_long_name = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_LONG_NAME")
out_dv_attr_var_desc = getenv("GFED_REGRID_HOURLY_EMIS_DATAVAR_ATTR_VAR_DESC")

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

;;; Create attributes for output data variable.
;;; Values below copy/mod from FIXME.
undef("get_output_var_attributes") ; allow redefinition in console
;  dv [*] [*] [*] [*] : float ; datavar
function get_output_var_attributes(
  dv : numeric ; datavar
)
begin
  dv@units = out_dv_attr_units
  dv@long_name = out_dv_attr_long_name
  dv@var_desc = out_dv_attr_var_desc
  return(dv)
end ; function get_output_var_attributes

;;; Create file/global attributes for output IOAPI.
;;; Values below copied mostly verbatim from FIXME.
undef("get_output_file_attributes") ; allow redefinition in console
function get_output_file_attributes()
local out_file_att, var_list_name
begin
  var_list_name = "VAR-LIST" ; IOAPI brainfart below
  out_file_att = True ; create attribute container
  out_file_att@IOAPI_VERSION = "$Id: @(#) ioapi library version 3.0 $                                           " ;
  out_file_att@EXEC_ID = "????????????????                                                                " ;
  out_file_att@FTYPE = 1 ;
  out_file_att@CDATE = 2011355 ;
  out_file_att@CTIME = 184611 ;
  out_file_att@WDATE = 2011355 ;
  out_file_att@WTIME = 184611 ;
  out_file_att@SDATE = 2008001 ;
  out_file_att@STIME = 0 ;
  out_file_att@TSTEP = 10000 ;
  out_file_att@NTHIK = 1 ;
  out_file_att@NCOLS = 459 ;
  out_file_att@NROWS = 299 ;
  out_file_att@NLAYS = 1 ;
  out_file_att@NVARS = 1 ; 52 in original
  out_file_att@GDTYP = 2 ;
  out_file_att@P_ALP = 33. ;
  out_file_att@P_BET = 45. ;
  out_file_att@P_GAM = -97. ;
  out_file_att@XCENT = -97. ;
  out_file_att@YCENT = 40. ;
  out_file_att@XORIG = -2556000. ;
  out_file_att@YORIG = -1728000. ;
  out_file_att@XCELL = 12000. ;
  out_file_att@YCELL = 12000. ;
  out_file_att@VGTYP = -1 ;
  out_file_att@VGTOP = tofloat(0) ;
  out_file_att@VGLVLS = tofloat( (/ 0,0 /) )
  out_file_att@GDNAM = "12US1_459X299   " ;
  out_file_att@UPNAM = "MRGGRID         " ;
;  out_file_att@VAR-LIST = "N2O             "
;  out_file_att@"VAR-LIST" = "N2O             "
  ; why did IOAPI put a dash in there?
  out_file_att@$var_list_name$ = "N2O             "
; out_file_att@FILEDESC = ; loonnngggg string, omitted
  out_file_att@HISTORY = out_attr_file_history
  return(out_file_att)
end ; function get_output_file_attributes

;;; Converts 3-hourly fractions (dims=(8, row, col))
;;; into hourly fractions (dims=(24, row, col))
undef("get_hourly_frac_arr") ; allow redefinition in console
function get_hourly_frac_arr(
  hour3ly_frac_arr [*] [*] [*] : numeric ; ASSERT: constant! won't be changed
)
local hourly_frac_arr, hourly_frac_arr_dims
begin
  hourly_frac_arr_dims = dimsizes(hour3ly_frac_arr)
  hourly_frac_arr_dims(0) = hours_per_day
  hourly_frac_arr = new (hourly_frac_arr_dims, "float") ; for the month
  do i_hour = 0, toint(hours_per_day - 1)
    reference_hour = toint(floor(i_hour / 3)) ; 3-hour period to use
    ; linearly interpolate 3-hour intervals -> 1-hour intervals
    hourly_frac_arr(i_hour,:,:) = hour3ly_frac_arr(reference_hour,:,:) / 3.0
  end do ; iterating hours
  ; convert units from .../hr to .../s
  hourly_frac_arr = hourly_frac_arr / seconds_per_hour
  return(hourly_frac_arr)
end ; function get_hourly_frac_arr

;;; Get the (string) path to the file for the hourly emissions for the day.
undef("get_daily_emis_fp") ; allow redefinition in console
function get_daily_emis_fp(
  month [1] : numeric, ; int index of month in year (1..12)
  day [1] : numeric,   ; int index of day in month (1..31)
  fp_template [1] : string,  ; template for file path
  date_template [1] : string ; template for date portion of file path
)
local date_string, fp
begin
;   date_string = sprinti("%0.2i%0.2i", month, day) ; '%0.'==left pad zeros
  ; NCL sprint* are so weak
  date_string = sprinti("%0.2i", month)
  date_string = date_string + sprinti("%0.2i", day)
  fp = str_sub_str(fp_template, date_template, date_string)
  return(fp)
end ; function get_daily_emis_fp

;;; Create output datavar with dimensions similar to those of a template,
;;; excepting the first 2. OK, so it's not so general :-)
undef("get_output_datavar") ; allow redefinition in console
function get_output_datavar(
  template_dv [*] [*] [*] [*] : numeric, ; template datavar and ...
  template_dv_name [1] : string,         ; ... its name, and ...
  template_dv_fh [1] : file              ; ... its file's handle
)
local dv, dv_dims, template_dv_dim_names, i_name
begin
  ;;; create raw output datavar
  dv_dims = dimsizes(template_dv)
  dv_dims(0) = out_tstep_dim_n    ; global
  dv_dims(1) = out_lay_dim_n      ; global
  dv = new(dv_dims, "float")

  ;;; transfer datavar coordvars
  ;; transfer dimension names
  template_dv_dim_names = \
    getfilevardims(template_dv_fh, template_dv_name) ; latter global
  do i_name = 0, 3 ; 4D
    dv!i_name = template_dv_dim_names(i_name)
  end do
  ;; transfer coordvar values? do I need this?
  dv&$template_dv_dim_names(0)$ = ispan(1,out_tstep_dim_n,1)
  dv&$template_dv_dim_names(1)$ = ispan(1,out_lay_dim_n,1)
  dv&$template_dv_dim_names(2)$ = template_dv&$template_dv_dim_names(2)$
  dv&$template_dv_dim_names(3)$ = template_dv&$template_dv_dim_names(3)$
  ;; transfer coordvar attributes? do I need this?
;   dv&$template_dv_dim_names(0)$ = \
;     get_tstep_coord_attributes(dv&$template_dv_dim_names(0)$)
;   dv&$template_dv_dim_names(1)$ = \
;     get_lay_coord_attributes(dv&$template_dv_dim_names(1)$)

  ;;; transfer datavar attributes
  dv = get_output_var_attributes(dv)
  return(dv)
end ; function get_output_datavar

;;; function definition template
; undef("%%%") ; allow redefinition in console
; function %%%(
;   arg [dim] : type, ; comment
; )
; local
; begin
;   return()
; end ; function %%%

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

begin ; skip if copy/paste-ing to console

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

  ;;; open regridded emissions and fractions
  ;; monthly emissions
  monthly_emis_fh = addfile(monthly_emis_fp, "r")
  monthly_emis_dv = monthly_emis_fh->$monthly_emis_dv_name$
;   monthly_emis_arr = (/ monthly_emis_dv /)
;   monthly_emis_dims = dimsizes(monthly_emis_arr) ; (/ 12, 299, 459 /)

  ;; daily fractions (by day for year)
  daily_frac_fh = addfile(daily_frac_fp, "r")
  daily_frac_dv = daily_frac_fh->$daily_frac_dv_name$
;  daily_frac_arr = (/ daily_frac_dv /)

  ;; 3-hourly fractions (by month for year)
  hour3ly_frac_fh = addfile(hour3ly_frac_fp, "r")
  hour3ly_frac_dv = hour3ly_frac_fh->$hour3ly_frac_dv_name$
;  hour3ly_frac_arr = (/ hour3ly_frac_dv /)

  ;;; setup daily emissions output
  ;; create output attributes
  out_file_atts = get_output_file_attributes()
;   ; start debugging---------------------------------------------------
;   printVarSummary(out_file_atts)
;   ;   end debugging---------------------------------------------------

  ;; Create output datavar with dimensions=(TSTEP | 25, LAY | 1, ROW | 299, COL | 459)
  ;; using 3hourly fractions as a template (since that has very similar schema)
  out_dv = get_output_datavar(\
    hour3ly_frac_dv, hour3ly_frac_dv_name, hour3ly_frac_fh)

;   ; start debugging---------------------------------------------------
;   printVarSummary(out_dv) ; check dimensions
; ; Variable: out_dv
; ; 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: [1..1]
; ;             ROW: [1854000..-1722000]
; ;             COL: [-2550000..2946000]
; ; Number Of Attributes: 4
; ;   _FillValue :  9.96921e+36
; ;   units :       moles/s         
; ;   long_name :   N2O             
; ;   var_desc :    N2O                                                                             
; ; TODO: why are ROW decreasing? and does it matter?
;   ;   end debugging---------------------------------------------------

  ;;; create unit-conversion matrix to finish mapping
  ;;; from GFED-supplied 3-hourly output (units==gN2O/m^2/(3-hour) per ftp://zea.ess.uci.edu/GFEDv3.1/Readme.pdf )
  ;;; to CMAQ-desired molN2O/s . Full conversion is

  ;;; gN2O        (3-hr)   m^2   molN2O   hr
  ;;; ----------* ------ * --- * ------ * --
  ;;; m^2 (3-hr)   3 hr          gN2O     s 

  ;;; but
  ;;; * I do the 3-hr-interval -> sec conversion in function get_hourly_frac_arr
  ;;; * I do the flux-rate -> molar-rate conversion in convert_monthlies.ncl
  ;;; so the following is unnecessary--wish I'd remembered that this morning :-(

;   out_areas_fh = addfile(out_areas_fp, "r")
;   out_areas_dv = out_areas_fh->$out_areas_dv_name$
;   out_conv_mx = (/ out_areas_dv /) / (molar_mass_N2O * seconds_per_hour)

;----------------------------------------------------------------------
; main loop(s): month, day, hour
;----------------------------------------------------------------------

  ;;; for each month in *_monthly_emis_emissions_regrid.nc
;  do n_month = 1, 2 ; for testing
;  do n_month = 12, 12 ; for debugging
  do n_month = 1, 12
    ; TODO: show progress
    i_month = n_month - 1 ; array pointer
    ;;; get the month's emissions
    monthly_emis_arr = (/ monthly_emis_dv(i_month,:,:) /)
    ; debugging: check dimensions == (/ 299, 459 /)
    ;;; get the month's 3hourly fractions
    hour3ly_frac_arr = (/ hour3ly_frac_dv(i_month,:,:,:) /)
    ; debugging: check dimensions == (/ 8, 299, 459 /)
    ;;; create "hourly fractions" array for the month
    hourly_frac_arr = get_hourly_frac_arr(hour3ly_frac_arr)

;     ; start debugging-------------------------------------------------
;     printVarSummary(hourly_frac_arr) ; check dimensions == (/ 24, 299, 459 /)
; ; Variable: hourly_frac_arr
; ; Type: float
; ; Total Size: 13175136 bytes
; ;             3293784 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes: [24] x [299] x [459]
; ; Coordinates: 
; ; Number Of Attributes: 1
; ;   _FillValue :        9.96921e+36
;     ;   end debugging-------------------------------------------------

    ;;; for each day in the month in *_daily_frac_fractions_regrid.nc
    ;;; (encoded as sequential int's, i.e., 1 Feb -> 32),
    ;;; write the hourly emissions for the day
    first_day_this_month = \ ; will need this later
      first_day_of_month(model_year, n_month)
;     last_day_this_month = \  ; debugging
;       last_day_of_month(model_year, n_month)
    do i_day = first_day_this_month, last_day_of_month(model_year, n_month)

      ;;; Setup daily emissions file path, or rather paths:
      ;;; to get CMAQ extension, gotta work around {CMAQ odd practice, NCL supported extensions}
      n_day = i_day - first_day_this_month + 1
      out_fp_cmaq = get_daily_emis_fp(\
        n_month, n_day, out_emis_fp_template_cmaq, day_template_string)
      out_fp_ncl = get_daily_emis_fp(\
        n_month, n_day, out_emis_fp_template_ncl, day_template_string)
      ;;; don't remake file if already present
      if      (isfilepresent(out_fp_cmaq)) then
        print(this_fn+": skipping output file='"+out_fp_cmaq+"'")
        continue ; skip this day
      else if (isfilepresent(out_fp_ncl)) then
        print(this_fn+": skipping output precursor='"+out_fp_ncl+"'")
        continue ; skip this day
      end if
      end if ; yes, again! see http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If

      ;; else, 
      ;; get daily fractions
      daily_frac_arr = (/ daily_frac_dv((i_day-1),:,:) /) ; C-style indexing!
;       ; start debugging-----------------------------------------------
;       printVarSummary(daily_frac_arr) ; check dimensions
; ; Variable: daily_frac_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 :        -3.4e+38
;       ;   end debugging-----------------------------------------------

      ;; compute daily emissions
      daily_emis_arr = monthly_emis_arr * daily_frac_arr
;       ; start debugging-----------------------------------------------
;       printVarSummary(daily_emis_arr) ; check dimensions
; ; Variable: daily_emis_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 :        -3.4e+38
;       ;   end debugging-----------------------------------------------

      ;; for each hour in the day, write hourly emissions to datavar
      do n_hour = 1, hours_per_day
        ; TODO: show progress
        i_hour = n_hour - 1
        ;; compute hourly emissions in molN2O/s : note conversions above
        hourly_emis_arr = daily_emis_arr * hourly_frac_arr(i_hour,:,:)
          
;         ; start debugging-----------------------------------------------
;         printVarSummary(hourly_emis_arr) ; check dimensions
; ; Variable: hourly_emis_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 :	-3.4e+38
;         ;   end debugging-----------------------------------------------

        ;; write daily emissions datavar (assuming |LAY|==1)
;        out_dv(i_hour,0,:,:) = hourly_emis_arr
; whacks {dimension names, coordinate variables} ROW, COL
        out_dv(i_hour,0,:,:) = (/ hourly_emis_arr /)
      end do ; iterating hours

      ;; since CMAQ wants one more TSTEP than actual hours, duplicate last hour
      out_dv(hours_per_day,0,:,:) = (/ out_dv(hours_per_day-1,0,:,:) /)

;       ; start debugging-----------------------------------------------
;       printVarSummary(out_dv) ; check dimensions
; ; Variable: out_dv
; ; 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: [1..1]
; ;             ROW: [1854000..-1722000]
; ;             COL: [-2550000..2946000]
; ; Number Of Attributes: 4
; ;   _FillValue :        9.96921e+36
; ;   units :     moles/s         
; ;   long_name : N2O             
; ;   var_desc :  N2O                                                                             
;       ;   end debugging-----------------------------------------------

      ;;; Write output daily emissions file to NCL filepath
      out_fh = addfile(out_fp_ncl, "c") ; create, not write
      ;; write global attributes: must do this before writing datavar
      out_file_atts@creation_date = systemfunc("date")
      fileattdef(out_fh, out_file_atts)
      ;; write datavar
      ; TODO: replace with real progress control!
      print("writing data to out_fp='"+out_fp_ncl+"'") 
      out_fh->$out_dv_name$ = out_dv
      ;; replace with desired CMAQ filepath
      system("mv "+out_fp_ncl+" "+out_fp_cmaq)

    end do ; iterating days
  end do ; iterating months

end ; make_hourlies.ncl