Source

MOZART-global to AQMEII-NA / hsp_to_gh_amsl__MOZART.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
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
;!/usr/bin/env ncl ; NCL::write_table requires version >= 6.1.0
;----------------------------------------------------------------------
; Convert MOZART hybrid sigma-pressure (hσp) grids (dimensions=(hσp, lat, lon))
; to the more "Cartesian" grids (at least vertically)

; double z_int_geo(ilev, lat, lon) # layer-interface elevations on geographical coordinates
; double z_mid_geo(lev, lat, lon)  # layer-midpoint elevations on geographical coordinates

; where z=geometric height above mean sea level.

; TODO: refactor this and hsp_to_gh_amsl__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 delete_VarAtts

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

; TODO: prototype args
function tdz_to_z(
  tdz_arr,     ; array returned by stdatmus_p2tdz
  level_name,  ; name for the level coordvar (e.g., "ilev", "lev")
  level_file,  ; source for vertical/level dimension
  lonlat_file  ; source for horizontal dimensions
)
  ; returns z_arr, an array with dimensions (level, lat, lon)
  local z_arr_dims, level_min, level_max, lat_min, lat_max, lon_min, lon_max
begin

  ; name the dimensions (which != coordinate variables!), use to collapse dimension=tdz
  tdz_arr!0 = "tdz"
  tdz_arr!1 = level_name ; a coordinate, not an integer index (like dim=lev)
  tdz_arr!2 = "lat"
  tdz_arr!3 = "lon"

  ; create elevation "variable" (to be)
;  z_arr = tdz_arr( tdz | 2, ilev | :, lat | :, lon | :)
  z_arr = tdz_arr( tdz | 2, $level_name$ | :, lat | :, lon | :)
  ; fix coordinates and attributes. TODO: copy from source grid(s)
  z_arr@_FillValue = todouble(tdz_arr@_FillValue)
  z_arr!0 = level_name
  z_arr!1 = "lat"
  z_arr!2 = "lon"
  z_arr&$level_name$ = level_file->$level_name$ ; not in sg_n2o_file
  z_arr&lat = lonlat_file->lat          ; but sg_other_file is not cropped
  z_arr&lon = lonlat_file->lon
  z_arr@long_name = "elevation at layer interfaces, computed by ncl::stdatmus_p2tdz"
  z_arr@units = "m"

  ; reverse level coord to make it increase with elevation: thanks DJS!
  z_arr = z_arr(::-1,:,:)
  z_arr_dims = dimsizes(z_arr)
  level_min = 0
  level_max = z_arr_dims(0)-1
  lat_min = 0
  lat_max = z_arr_dims(1)-1
  lon_min = 0
  lon_max = z_arr_dims(2)-1

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

; ; using layer midpoints, HEIGHT IS NOT A COORDINATE!
; ; we have NA heights ... and height varies by gridcell ...

;   printVarSummary(z_arr)
; ; Variable: z_arr
; ; Type: float
; ; Total Size: 143640 bytes
; ;             35910 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes: [ilev | 57] x [lat | 21] x [lon | 30]
; ; Coordinates:
; ;             ilev: [1000..1.650789985433221]
; ;             lat: [59.68421052631579..21.78947368421052]
; ;             lon: [-130..-57.5]
; ; Number Of Attributes: 3
; ;   units :     m
; ;   long_name : elevation at layer interfaces, computed by ncl::stdatmus_p2tdz
; ;   _FillValue :        -999

;   ; "do the corners" for z_arr: do [top,bot] for [SW,NW,NE,SE]

;   print("") ; can't print("\n")
;   print("z_arr(min("+level_name+"),min(lat),min(lon))==")
;   print(z_arr(level_min,lat_min,lon_min))
; ; (0)   1157.416076155514  ; for z_int

;   print("")
;   print("z_arr(max("+level_name+"),min(lat),min(lon))==")
;   print(z_arr(level_max,lat_min,lon_min))
; ; (0)   43899.81501286371  ; for z_int

;   print("")
;   print("z_arr(min("+level_name+"),max(lat),min(lon))==")
;   print(z_arr(level_min,lat_max,lon_min))
; ; (0)   -999               ; for z_int

;   print("")
;   print("z_arr(max("+level_name+"),max(lat),min(lon))==")
;   print(z_arr(level_max,lat_max,lon_min))
; ; (0)   43899.81501286371  ; for z_int

;   print("")
;   print("z_arr(min("+level_name+"),max(lat),max(lon))==")
;   print(z_arr(level_min,lat_max,lon_max))
; ; (0)   -999               ; for z_int

;   print("")
;   print("z_arr(max("+level_name+"),max(lat),max(lon))==")
;   print(z_arr(level_max,lat_max,lon_max))
; ; (0)   43899.81501286371  ; for z_int

;   print("")
;   print("z_arr(min("+level_name+"),min(lat),max(lon))==")
;   print(z_arr(level_min,lat_min,lon_max))
; ; (0)   5.618337654996648  ; for z_int

;   print("")
;   print("z_arr(max("+level_name+"),min(lat),max(lon))==")
;   print(z_arr(level_max,lat_min,lon_max))
; ; (0)   43899.81501286371  ; for z_int

; ;   end debug---------------------------------------------------------
  return(z_arr)
end ; function tdz_to_z

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

begin ; skip if copy/paste-ing to console
  ;; constants--------------------------------------------------------

  ; source grids

  ; MOZART-derived file with source datavar = N2O
  sg_n2o_fp = getenv("SG_N2O_FP")
  sg_n2o_file = addfile(sg_n2o_fp,"r")

  ; MOZART-derived file with source datavar = PS
  sg_ps_fp = getenv("SG_PS_FP")
  sg_ps_file = addfile(sg_ps_fp,"r")

  ; MOZART-derived file with all other datavars
  sg_other_fn = getenv("SG_OTHER_FN")
  sg_other_fp = getenv("SG_OTHER_FP")
  sg_other_file = addfile(sg_other_fp,"r")

  ; ETOPO1-derived file of surface elevations (only)
  sg_etopo1_fp = getenv("SG_ETOPO1_FP")
  sg_etopo1_file = addfile(sg_etopo1_fp,"r")

  ; target grids
  datavar_z_int_name = getenv("DATAVAR_Z_INT_NAME") ; datavar for elevations of layer interfaces
  datavar_z_mid_name = getenv("DATAVAR_Z_MID_NAME") ; datavar for elevations of layer midpoints

  ; target ASCII table file: only z_int (for data interchange with R)
  tg_ascii_fp = getenv("TG_ASCII_FP")
  ; will get handle for writing later
  tg_ascii_na = getenv("TG_ASCII_NA")
;  tg_ascii_line_format = getenv("TG_ASCII_LINE_FORMAT")
  tg_ascii_datum_format = getenv("TG_ASCII_DATUM_FORMAT")

  ; target netCDF: includes all datavars above, plus z_int_geo and z_mid_geo
  tg_netcdf_fp = getenv("TG_NETCDF_FP")
  ; will get handle for writing later

  ; coordinate dimensions and indices (used later for iteration)
  ilev_var = sg_other_file->ilev ; not in sg_n2o_file
  ilev_min = 0                ; C-style, min array index
  ilev_n = dimsizes(ilev_var)
  ilev_max = ilev_n - 1       ; max array index

  lev_var = sg_other_file->lev
  lev_min = 0
  lev_n = dimsizes(lev_var)
  lev_max = lev_n - 1

  lon_var = sg_n2o_file->lon     ; sg_other_file is not cropped (but probably does not matter)
  lon_min = 0
  lon_n = dimsizes(lon_var)
  lon_max = lon_n - 1

  lat_var = sg_n2o_file->lat     ; sg_other_file is not cropped (but probably does not matter)
  lat_min = 0
  lat_n = dimsizes(lat_var)
  lat_max = lat_n - 1

  ; magic
  ; index of z in dim=1 of stdatmus_p2tdz(p)
  tdz_height_i = 2
  ; stdatmus_p2tdz internally uses NA==-999 ???
  stdatmus_p2tdz_NA = -999

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

  N2O_var = sg_n2o_file->N2O ; draws
;> warning:FileReadVarAtt: _FillValue attribute type differs from variable (N2O) type in file (2008N2O_restart_region_N2O), forcing type conversion; may result in overflow and/or loss of precision
  ; try to fix warning
  ; TODO: don't hardcode type!
  N2O_var@_FillValue = todouble(N2O_var@_FillValue)
  N2O_var@missing_value = todouble(N2O_var@_FillValue)

  PS_var = sg_ps_file->PS ; draws
  ; try to fix warning
  PS_var@_FillValue = todouble(PS_var@_FillValue)
  PS_var@missing_value = todouble(PS_var@_FillValue)

  hyam_var = sg_other_file->hyam
  hybm_var = sg_other_file->hybm
  ; instead, do interfaces (see problem below with layer midpoint heights)
  hyai_var = sg_other_file->hyai
  hybi_var = sg_other_file->hybi

  ; for completeness, written below
  P0_var = sg_other_file->P0
  date_var = sg_other_file->date
  datesec_var = sg_other_file->datesec

  ; convert hσp -> pressure-------------------------------------------

  ; hσp -> pressure: "mid"=layer midpoints, "int"=layer interfaces
  P_int_arr = pres_hybrid_ccm(PS_var, P0_var, hyai_var, hybi_var)
  P_mid_arr = pres_hybrid_ccm(PS_var, P0_var, hyam_var, hybm_var)

  ; stdatmus_p2tdz wants mb! though it says nothing on its webpage :-(
  P_int_arr_in_mb = P_int_arr/100.0 ; Pa -> hPa
  P_mid_arr_in_mb = P_mid_arr/100.0

  ; for correcting stdatmus_p2tdz NAs (below)
  zS_arr = sg_etopo1_file->zS
  ; try to fix warning
  zS_arr@_FillValue = tofloat(zS_arr@_FillValue)
  zS_arr@missing_value = tofloat(zS_arr@_FillValue)

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

;   printVarSummary(N2O_var)                          
; ; Variable: N2O_var
; ; Type: float
; ; Total Size: 141120 bytes
; ;             35280 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes: [lev | 56] x [lat | 21] x [lon | 30]
; ; Coordinates:
; ;             lev: [1..56]
; ;             lat: [59.68421052631579..21.78947368421052]
; ;             lon: [-130..-57.5]
; ; Number Of Attributes: 8
; ;   _FillValue :        -1
; ;   units :     ppmV
; ;   missing_value :       -1
; ;   long_name : N2O concentration
; ;   projection :        +proj=longlat +datum=WGS84
; ;   projection_format : PROJ.4
; ;   min :       <ARRAY of 56 elements>
; ;   max :       <ARRAY of 56 elements>

;   printVarSummary(P_int_arr_in_mb)
; ; Variable: P_int_arr_in_mb
; ; Type: float
; ; Total Size: 143640 bytes
; ;             35910 values
; ; Number of Dimensions: 3
; ; Dimensions and sizes: [57] x [21] x [30]
; ; Coordinates:
; ; Number Of Attributes: 1
; ;   _FillValue :        -1

;   ; do [top,bot] for [SW,NW,NE,SE]
;   print(P_int_arr_in_mb(ilev_min,lat_min,lon_min))
; ; (0)   1.65079
;   print(P_int_arr_in_mb(ilev_max,lat_min,lon_min))
; ; (0)   881.7172
;   print(P_int_arr_in_mb(ilev_min,lat_max,lon_min))
; ; (0)   1.65079
;   print(P_int_arr_in_mb(ilev_max,lat_max,lon_min))
; ; (0)   1016.549
;   print(P_int_arr_in_mb(ilev_min,lat_max,lon_max))
; ; (0)   1.65079
;   print(P_int_arr_in_mb(ilev_max,lat_max,lon_max))
; ; (0)   1021.245
;   print(P_int_arr_in_mb(ilev_min,lat_min,lon_max))
; ; (0)   1.65079
;   print(P_int_arr_in_mb(ilev_max,lat_min,lon_max))
; ; (0)   1012.575

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

  ; pressure -> geometric height AMSL (hopefully)
  tdz_int_arr = stdatmus_p2tdz(P_int_arr_in_mb) ;
  tdz_int_arr@_FillValue = todouble(stdatmus_p2tdz_NA)
  tdz_int_arr@missing_value = todouble(stdatmus_p2tdz_NA)

  tdz_mid_arr = stdatmus_p2tdz(P_mid_arr_in_mb) ;
  tdz_mid_arr@_FillValue = todouble(stdatmus_p2tdz_NA)
  tdz_mid_arr@missing_value = todouble(stdatmus_p2tdz_NA)

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

;  print("") ; newline
;  print("printVarSummary(tdz_int_arr)==")
;  printVarSummary(tdz_int_arr)
;; Variable: tdz_int_arr
;; Type: float
;; Total Size: 430920 bytes
;;             107730 values
;; Number of Dimensions: 4
;; Dimensions and sizes: [3] x [57] x [21] x [30]
;; Coordinates:
;; Number Of Attributes: 1
;;   _FillValue :        -999
;
;  print("") ; newline
;  print("printVarSummary(tdz_mid_arr)==")
;  printVarSummary(tdz_mid_arr)

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

  ; create elevation datavars from stdatmus_p2tdz output
  z_int_geo_var = tdz_to_z(tdz_int_arr, "ilev", sg_other_file, sg_n2o_file)
  z_mid_geo_var = tdz_to_z(tdz_mid_arr, "lev", sg_other_file, sg_n2o_file)
  ; but we'll compute midpoints later, instead ...

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

  print("") ; newline
  print("printVarSummary(z_int_geo_var)==")
  printVarSummary(z_int_geo_var)

  print("") ; newline
  print("printVarSummary(z_mid_geo_var)==")
  printVarSummary(z_mid_geo_var)

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

  ; fix z_int_geo_var values------------------------------------------
  ; replace NA/_FillValue with ETOPO1 data if @ surface/min(ilev)

  ; find all locations with z=NA: replace if surface, throw if not surface
  m = 0 ; their number
  n = 0 ; total obs in array (yes, I know ...)
  do i=ilev_min, ilev_max
    do j=lat_min, lat_max
      do k=lon_min, lon_max
        n = n + 1
        ; see warning regarding NA @ http://www.ncl.ucar.edu/Document/Manuals/Ref_Manual/NclStatements.shtml#If
        if (ismissing(z_int_geo_var(i,j,k))) then
          m = m + 1
          ; are all NA locations ilev=0, i.e., surface?
          if (i .eq. ilev_min) then
            ; replace with corresponding value from ETOPO1 regrid,
            ; and add that value to all vertical layers at that horizontal location
            rep = zS_arr(j,k)        ; same (lat, lon)
            if (ismissing(rep)) then
              ; start debugging---------------------------------------
              print("ERROR: both z_int_geo_var and zS_arr NA @ ("+i+","+j+","+k+")")
              ;   end debugging---------------------------------------
              ; <sigh/> above is Caribbean, so replace with 0
              rep = 0
            end if ; (ismissing(rep))
            ; don't replace missing value with benthic/negative elevations
            ; TODO: "correct" ETOPO1 s.t. oceanic elevations=0 (sea level)
            ; since following approach is error for, e.g., Death Valley
            if (rep .lt. 0) then
              rep = 0
            end if

            ; stdatmus_p2tdz outputs an altitude, so one must
            ; add the surface elevation to each layer, not just the first.
;             ; start debugging---------------------------------------
;             print("WARNING: adding "+rep+" to all z(:,"+j+","+k+")")
;             ;   end debugging---------------------------------------
            z_int_geo_var(0,j,k) = 0
            z_int_geo_var(:,j,k) = rep + z_int_geo_var(:,j,k)
          else ; if !(i .eq. ilev_min)
            print("ERROR: non-surface elevation NA @ ("+i+","+j+","+k+")")
            ; fortunately none of those! otherwise, we're hosed
          end if ; (i .eq. ilev_min), i.e., surface
        end if ; (ismissing(z_int_geo_var(i,j,k)))

      end do ; k=lon_min, lon_max
    end do ; j=lat_min, lat_max
  end do ; i=ilev_min, ilev_max
  print("|NA in z_int_geo_var|="+m+", |z_int_geo_var|="+n)
; (0)   |NA in z_int_geo_var|=233, |z_int_geo_var|=35910

  ; write z_int_geo_var to ASCII table--------------------------------

  ; wish I could use
  ; write_table(tg_ascii_fp, "w",  ; overwrite tg_ascii_fp
  ;   ; list of values to write
  ;   ; format for writing values
  ; )
  ; but don't see how to go from z(lat,lon) to NCL-style list. So instead:
  ; one more pass, over z_int_geo_var: I'll take correctness over performance, for now

  ; vector=lines holds strings of form [lon, lat, z_int]
  lines_i = 0 ; index into vector=lines
  ; lines_n == |records| in vector=lines == length(z_int_geo_var)
  lines_n = ilev_n * lat_n * lon_n
  lines_max = lines_n - 1
  lines = new( (/ lines_n /), string)
  lines@_FillValue = tg_ascii_na

  ; loop order (hopefully) facilitates reading elevations vector @ each horizontal
  do j=lat_min, lat_max
    do k=lon_min, lon_max
      do i=ilev_min, ilev_max
;       lines(lines_i) = sprintf(tg_ascii_line_format, lon_var(k), lat_var(j), z_int_geo_var(i,j,k))
;       lines(lines_i) = sprintf(tg_ascii_line_format, (/ lon_var(k), lat_var(j), z_int_geo_var(i,j,k) /) )
        ; that would be too easy !-(
	line = sprintf(tg_ascii_datum_format, lon_var(k)) + "," + \
               sprintf(tg_ascii_datum_format, lat_var(j)) + "," + \
               sprintf(tg_ascii_datum_format, z_int_geo_var(i,j,k))
        lines(lines_i) = line
        lines_i = lines_i + 1
      end do ; i=ilev_min, ilev_max
    end do ; k=lon_min, lon_max
  end do ; j=lat_min, lat_max

;; start debugging---------------------------------------
;;  print("length(lines)=="+lines_i+" (should=="+lines_n+")")
;  printVarSummary(lines)
;  print("head(lines)==")
;  print(lines(ilev_min:ilev_max))
;  print("tail(lines)==")
;;  print(lines(lines_max))
;  print(lines((lines_max-ilev_max):lines_max))
;;   end debugging---------------------------------------

  asciiwrite (tg_ascii_fp, lines)

;; start debugging---------------------------------------
;  print("head("+tg_ascii_fp+")==")
;  system("head "+tg_ascii_fp) 
;  print("tail("+tg_ascii_fp+")==")
;  system("tail "+tg_ascii_fp) 
;;   end debugging---------------------------------------

  ; populate z_mid_geo_var values-------------------------------------
  ; assume z(midpoint(layer(i)))==mean(z_int(layer(i)), z_int(layer(i+1)))
  ; second pass is inefficient but avoids complications

  do j=lat_min, lat_max
    do k=lon_min, lon_max
      ilevs = z_int_geo_var(:,j,k) ; vertical "slice" @ horizontal location
      do i=lev_min, lev_max ; not ilev
        lo = ilevs(i)   ; elevation immediately below midpoint
        hi = ilevs(i+1) ; elevation immediately above midpoint
        if      (ismissing(lo)) then
          ; TODO: throw error and exit
          print("ERROR: NA @ z_int_geo_var("+i+","+j+","+k+")")
        else if (ismissing(hi)) then
          print("ERROR: NA @ z_int_geo_var("+(i+1)+","+j+","+k+")")
        else if (hi .le. lo) then
          print("ERROR: z_int not monotonic @ ("+i+","+j+","+k+")")
          ; start debugging---------------------------------------
          print("  z_int("+i+","+j+","+k+")=="+lo)
          print("  z_int("+(i+1)+","+j+","+k+")=="+hi)
          ;   end debugging---------------------------------------
        else
          z_mid_geo_var(i,j,k) = avg( (/ lo, hi /) )
        end if ; note NCL's goofy "'else if'==nested if"
        end if
        end if
      end do ; i=lev_min, lev_max
    end do ; k=lon_min, lon_max
  end do ; j=lat_min, lat_max

  ; write elevations (both) to netCDF file----------------------------
  ; using efficient/definition mode: see http://www.ncl.ucar.edu/Applications/method_2.shtml

  ; this will fail if tg_netcdf_fp exists!
  tg_file_z = addfile(tg_netcdf_fp, "c")
  setfileoption(tg_file_z, "DefineMode", True) ; enter define mode

  ; set global attributes
  tg_file_z_att               = True ; assign file attributes
  tg_file_z_att@title         = "layer interface elevations derived from "+sg_other_fn
  tg_file_z_att@source_file   =  sg_other_fn
  tg_file_z_att@Conventions   = "None"  
  tg_file_z_att@history = "see https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na"
  tg_file_z_att@creation_date = systemfunc("date")       
  fileattdef(tg_file_z, tg_file_z_att) ; copy file attributes   

  ; set global dimensions
  dimNames = (/ "lev", "lat", "lon", "ilev" /) ; same order as `ncdump -h`
  dimSizes = (/ lev_n, lat_n, lon_n, ilev_n /)
  dimUnlim = (/ False, False, False, False /)  ; no unlimited dimensions
  filedimdef(tg_file_z, dimNames, dimSizes, dimUnlim)

  ; define variables
  filevardef(tg_file_z, "lev", typeof(lev_var), getvardims(lev_var))
  filevardef(tg_file_z, "lat", typeof(lat_var), getvardims(lat_var))
  filevardef(tg_file_z, "lon", typeof(lon_var), getvardims(lon_var))   
  filevardef(tg_file_z, "ilev", typeof(ilev_var), getvardims(ilev_var))
  filevardef(tg_file_z, "N2O", typeof(N2O_var), getvardims(N2O_var))
  filevardef(tg_file_z, "PS", typeof(PS_var), getvardims(PS_var))
  filevardef(tg_file_z, "date", typeof(date_var), getvardims(date_var))
  filevardef(tg_file_z, "datesec", typeof(datesec_var), getvardims(datesec_var))
  filevardef(tg_file_z, "hyai", typeof(hyai_var), getvardims(hyai_var))
  filevardef(tg_file_z, "hybi", typeof(hybi_var), getvardims(hybi_var))
  filevardef(tg_file_z, "hyam", typeof(hyam_var), getvardims(hyam_var))
  filevardef(tg_file_z, "hybm", typeof(hybm_var), getvardims(hybm_var))
  filevardef(tg_file_z, "P0", typeof(P0_var), getvardims(P0_var))

;;; TODO: do function for both datavar={datavar_z_int_name (for int), datavar_z_mid_name (for mid)}
  filevardef(tg_file_z, datavar_z_int_name, typeof(z_int_geo_var), getvardims(z_int_geo_var))
  filevardef(tg_file_z, datavar_z_mid_name, typeof(z_mid_geo_var), getvardims(z_mid_geo_var))

  ; define variables' attributes
  ; copy all attributes for coordvars
  filevarattdef(tg_file_z, "lev", lev_var)
  filevarattdef(tg_file_z, "lat", lat_var)
  filevarattdef(tg_file_z, "lon", lon_var)
  filevarattdef(tg_file_z, "ilev", ilev_var)

  ; copy most attributes for horizontal datavars
  delete_VarAtts(N2O_var, (/ "min", "max" /))   ; may have changed
  filevarattdef(tg_file_z, "N2O", N2O_var)

  delete_VarAtts(PS_var, (/ "min", "max" /))    ; may have changed
  filevarattdef(tg_file_z, "PS", PS_var)

;   delete_VarAtts(z_int_geo_var, (/ "min", "max" /)) ; may have changed
;   delete_VarAtts(z_int_geo_var, (/ "lat", "lon" /)) ; why are these here?
  z_int_geo_var@long_name = "elevation at layer interfaces"
  filevarattdef(tg_file_z, datavar_z_int_name, z_int_geo_var)

;   delete_VarAtts(z_mid_geo_var, (/ "min", "max", "lat", "lon" /))
  z_mid_geo_var@long_name = "elevation at layer midpoints"
  filevarattdef(tg_file_z, datavar_z_mid_name, z_mid_geo_var)

  ; copy non-horizontal datavar attributes completely
  filevarattdef(tg_file_z, "date", date_var)
  filevarattdef(tg_file_z, "datesec", datesec_var)
  filevarattdef(tg_file_z, "hyai", hyai_var)
  filevarattdef(tg_file_z, "hybi", hybi_var)
  filevarattdef(tg_file_z, "hyam", hyam_var)
  filevarattdef(tg_file_z, "hybm", hybm_var)
  filevarattdef(tg_file_z, "P0", P0_var)

  setfileoption(tg_file_z, "DefineMode", False) ; exit define mode
  ; write data to defined file
  tg_file_z->N2O = (/ N2O_var /)
  tg_file_z->PS = (/ PS_var /)
  tg_file_z->$datavar_z_int_name$ = (/ z_int_geo_var /)
  tg_file_z->$datavar_z_mid_name$ = (/ z_mid_geo_var /)
  tg_file_z->date = (/ date_var /)
  tg_file_z->datesec = (/ datesec_var /)
  tg_file_z->hyai = (/ hyai_var /)
  tg_file_z->hybi = (/ hybi_var /)
  tg_file_z->hyam = (/ hyam_var /)
  tg_file_z->hybm = (/ hybm_var /)
  tg_file_z->P0 = (/ P0_var /)

end