1. Tom Roche
  2. MOZART-global to AQMEII-NA

Commits

Tom Roche  committed d55823a

Begin separating code, README into MOZART- and CMAQ-sides.

Eventually will need to refactor code, but for now I don't fully understand the similarities and differences of MOZART- and CMAQ-side function, so will just copy/mod: "it's the simplest thing that could possibly work."

TODO:
next major:
* plan item=3.3: read CSV back to NCL, convert (hσp,lat,lon) to (x,y,z).
minor:
* debug why writing CSV from NCL is so slow: netCDF datavar (probably most of 13 MB file) takes 19.22 min to write a 66 MB CSV
* debug why hsp_to_gh_amsl.ncl converts type(z_int) float->double
* debug why regrid.surface.elevation.r input.max(z) > output.max(z), when just regridding over same region
* ping r-sig-geo regarding raster datatype truncation
* refactor the now-6 sets of {.sh, .ncl, .r}
* investigate fringe of apparently-missing data along north and east borders of surface-elevation plot

  • Participants
  • Parent commits 9f9eca9
  • Branches master

Comments (0)

Files changed (5)

File README.md

View file
 
 1. Clone this repo.
 2. `cd` to its working directory.
-3. Visualize the global data:
-    1. Edit `./levelplot_netCDF_global.sh` to match your local configuration (e.g., set your PDF viewer).
-    2. Run `./levelplot_netCDF_global.sh` to visualize the raw data.
-    3. Check your resulting visualization against [the reference (download)][levelplot_netCDF_global.pdf].
-1. Create and visualize regional subsets of the global data:
-    1. Edit `./restrict_lon_lat_to_region.sh` to match your local configuration.
-    2. Run `./restrict_lon_lat_to_region.sh` to:
-        1. crop the global raw input data to the region (and save separately).
-        2. visualize the regional data.
-    3. Check your resulting visualizations against provided references (downloads):
-        * [`2008N2O_restart_region_N2O.pdf`][2008N2O_restart_region_N2O.pdf]
-        * [`2008N2O_restart_region_PS.pdf`][2008N2O_restart_region_PS.pdf]
-        * [`ETOPO1_Ice_g_gmt4.grd_region_zS.pdf`][ETOPO1_Ice_g_gmt4.grd_region_zS.pdf]
-1. Regrid the surface elevation data to match the grid of the other regional inputs:
-    1. Edit `./regrid_surface_elevation.sh` to match your local configuration.
-    2. Run `./regrid_surface_elevation.sh` to:
-        1. regrid the surface elevation data (and save separately).
-        2. visualize the regridded data.
-    3. Check your resulting visualizations against [the reference (download)][ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.pdf].
-1. Compute layer-interface and layer-midpoint elevations:
-    1. Edit `./hsp_to_gh_amsl.sh` to match your local configuration.
-    1. Run `./hsp_to_gh_amsl.sh`.
-    1. Check the values of the resulting new elevation datavars (e.g.,
-
-            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
-
-     ) in the resulting netCDF output against [the reference (download)][2008N2O_restart_region_z.nc]. (Yes, this needs visualizations ...)
-    1. Check the values of `z_int_geo` written (as ASCII CSV, for data interchange with R) to [the reference (download)][2008N2O_restart_region_z.txt]. (Yes, this also needs visualizations ...)
-
-1. Compute fully-Cartesian layer-interface elevations:
-    1. Edit `./cartesianize_elevations.sh` to match your local configuration.
-    1. Run `./cartesianize_elevations.sh`. Tested on both tlrPanP5 and terrae.
-    1. TODO: provide visualization to check the values of the resulting elevation datavar==`nodeCoords`
+3. Process MOZART-side data:
+    1. Visualize the global data:
+        1. Edit `./levelplot_netCDF_global.sh` to match your local configuration (e.g., set your PDF viewer).
+        2. Run `./levelplot_netCDF_global.sh` to visualize the raw data.
+        3. Check your resulting visualization against [the reference (download)][levelplot_netCDF_global.pdf].
+    1. Create and visualize regional subsets of the global data:
+        1. Edit `./restrict_lon_lat_to_region.sh` to match your local configuration.
+        2. Run `./restrict_lon_lat_to_region.sh` to:
+            1. crop the global raw input data to the region (and save separately).
+            2. visualize the regional data.
+        3. Check your resulting visualizations against provided references (downloads):
+            * [`2008N2O_restart_region_N2O.pdf`][2008N2O_restart_region_N2O.pdf]
+            * [`2008N2O_restart_region_PS.pdf`][2008N2O_restart_region_PS.pdf]
+            * [`ETOPO1_Ice_g_gmt4.grd_region_zS.pdf`][ETOPO1_Ice_g_gmt4.grd_region_zS.pdf]
+    1. Regrid the surface elevation data to match the grid of the other regional inputs:
+        1. Edit `./regrid_surface_elevation.sh` to match your local configuration.
+        2. Run `./regrid_surface_elevation.sh` to:
+            1. regrid the surface elevation data (and save separately).
+            2. visualize the regridded data.
+        3. Check your resulting visualizations against [the reference (download)][ETOPO1_Ice_g_gmt4.grd_region_zS_regrid.pdf].
+    1. Compute layer-interface and layer-midpoint elevations:
+        1. Edit `./hsp_to_gh_amsl__MOZART.sh` to match your local configuration.
+        1. Run `./hsp_to_gh_amsl__MOZART.sh`.
+        1. Check the values of the resulting new elevation datavars (e.g.,
+
+                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
+
+         ) in the resulting netCDF output against [the reference (download)][2008N2O_restart_region_z.nc]. (Yes, this needs visualizations ...)
+        1. Check the values of `z_int_geo` written (as ASCII CSV, for data interchange with R) to [the reference (download)][2008N2O_restart_region_z.txt]. (Yes, this also needs visualizations ...)
+
+    1. Compute fully-Cartesian layer-interface elevations:
+        1. Edit `./cartesianize_elevations.sh` to match your local configuration.
+        1. Run `./cartesianize_elevations.sh`. Tested on both tlrPanP5 and terrae.
+        1. TODO: provide visualization to check the values of the resulting elevation datavar==`nodeCoords`
+4. Process CMAQ-side data:
 
 [levelplot_netCDF_global.pdf]: https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/levelplot_netCDF_global.pdf
 [2008N2O_restart_region_N2O.pdf]: https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_N2O.pdf

File hsp_to_gh_amsl.ncl

-;!/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.
-;----------------------------------------------------------------------
-
-;----------------------------------------------------------------------
-; 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 ; hsp_to_gh_amsl.ncl

File hsp_to_gh_amsl.sh

-#!/usr/bin/env bash
-# Requires NCL.
-# Driver for hsp_to_gh_amsl.ncl--configure as needed for your platform.
-# Convert MOZART hybrid sigma-pressure grid (dimensions=(hsp, lat, lon))
-# to a more "Cartesian" grid with z=geometric height above mean sea level
-# (dimensions=(gh_AMSL, lat, lon)). Should work with older NCL versions.
-
-# constants with some simple manipulations------------------------------
-
-THIS_FN="$(basename $0)"
-# will do the real work
-
-# TODO: create from THIS_FN, e.g., use `sed -e s/_/./g`
-export CALL_NCL_FN='hsp_to_gh_amsl.ncl'
-
-# modules on EMVL. TODO: fixme! (these fail from here, are currently just reminders.)
-# You may not need them at all!
-module add netcdf-4.1.2
-#module add ncl_ncarg-6.0.0
-# Temporarily using version on terrae installed by Jerry Herwehe (thanks!)
-#module add nco-4.0.5
-
-# Ensure $NCARG_ROOT is set for NCL--it should be in your dotfiles, but, JIC
-# tlrPanP5
-#NCL_VERSION='ncl_ncarg-6.1.0.Linux_Debian_x86_64_nodap_gcc445'
-#export NCARG_ROOT="${HOME}/bin/${NCL_VERSION}"
-## Only change default if absolutely necessary!
-#NCL_EXEC='ncl'
-# terrae temporary (they gotta upgrade RHEL)
-NCL_VERSION='ncl_ncarg/ncl_ncarg-6.1.0.Linux_RedHat_x86_64_gcc412'
-export NCARG_ROOT="/home/hhg/${NCL_VERSION}"
-NCL_EXEC="${NCARG_ROOT}/bin/ncl"
-
-# workspace
-export WORK_DIR="$(pwd)" # keep it simple for now: same dir as top of repo
-# will do the real work (restricting and plotting)
-export CALL_NCL_FP="${WORK_DIR}/${CALL_NCL_FN}"
-
-# # for plotting
-# PDF_VIEWER='evince'  # whatever works on your platform
-# # temporally disaggregate multiple plots
-# DATE_FORMAT='%Y%m%d_%H%M'
-# export PDF_DIR="${WORK_DIR}"
-
-# source grids
-
-# directories
-SG_DIR="${WORK_DIR}" # keeping it simple
-TG_DIR="${WORK_DIR}"
-
-# files and URIs (from bitbucket repo)
-SG_ORIGINAL_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart.nc'                       # uncropped
-SG_ORIGINAL_FN="$(basename ${SG_ORIGINAL_URI})"
-SG_ORIGINAL_FN_ROOT="${SG_ORIGINAL_FN%.*}"
-SG_ORIGINAL_FN_EXT="${SG_ORIGINAL_FN##*.}"
-NETCDF_EXT="${SG_ORIGINAL_FN_EXT}"                 # stay consistent
-ASCII_EXT='txt'
-
-SG_CROPPED_FN_ROOT="${SG_ORIGINAL_FN_ROOT}_region" # cropped
-SG_CROPPED_FN="${SG_CROPPED_FN_ROOT}.${NETCDF_EXT}"
-
-# MOZART-derived file with source datavar=N2O
-SG_N2O_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_N2O.nc'
-SG_N2O_FN="$(basename ${SG_N2O_URI})"
-export SG_N2O_FP="${SG_DIR}/${SG_N2O_FN}"
-
-# MOZART-derived file with source datavar=PS
-SG_PS_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
-SG_PS_FN="$(basename ${SG_PS_URI})"
-export SG_PS_FP="${SG_DIR}/${SG_PS_FN}"
-
-# ETOPO1-derived file of surface elevations (only)
-SG_ETOPO1_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
-SG_ETOPO1_FN="$(basename ${SG_ETOPO1_URI})"
-export SG_ETOPO1_FP="${SG_DIR}/${SG_ETOPO1_FN}"
-
-# MOZART-derived file with all other datavars
-SG_OTHER_URI="${SG_ORIGINAL_URI}"
-export SG_OTHER_FN="$(basename ${SG_OTHER_URI})" # for doc strings
-export SG_OTHER_FP="${SG_DIR}/${SG_OTHER_FN}"
-
-# target grids
-
-# datavar for elevations of layer interfaces, but horizontally geographical
-export DATAVAR_Z_INT_NAME='z_int_geo'
-# datavar for elevations of layer midpoints, but [ditto]
-export DATAVAR_Z_MID_NAME='z_mid_geo'
-# TODO: (RSN) create fully Cartesian, horizontally Cartesian datavars ...
-
-# target netCDF file ~= ${SG_OTHER_FN} + new datavars above
-TG_NETCDF_FN="${SG_CROPPED_FN_ROOT}_z.${NETCDF_EXT}"
-export TG_NETCDF_FP="${TG_DIR}/${TG_NETCDF_FN}"
-
-# target ASCII file: write ${DATAVAR_Z_INT_NAME} for data interchange with R
-TG_ASCII_FN="${SG_CROPPED_FN_ROOT}_z.${ASCII_EXT}"
-export TG_ASCII_FP="${TG_DIR}/${TG_ASCII_FN}"
-# missing value for string output to be read by R
-export TG_ASCII_NA='NA,NA,NA'
-# not actually usable by @#$%^&! NCL::sprintf
-# export TG_ASCII_LINE_FORMAT='%f,%f,%f' # sprintf string for CSV line
-export TG_ASCII_DATUM_FORMAT='%f' # sprintf string for each datum in line
-
-# payload---------------------------------------------------------------
-
-# setup-----------------------------------------------------------------
-
-mkdir -p ${WORK_DIR}
-
-# get (or remove) data------------------------------------------------
-
-# NCL will not overwrite target ASCII, so we will. (TODO: backup/rename)
-if [[ -r "${TG_ASCII_FP}" ]] ; then
-  echo -e "WARNING: ${THIS_FN}: about to delete previously output ASCII file='${TG_ASCII_FP}'"
-  for CMD in \
-    "rm ${TG_ASCII_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# NCL will not overwrite target netCDF, so we will. (TODO: backup/rename)
-if [[ -r "${TG_NETCDF_FP}" ]] ; then
-  echo -e "WARNING: ${THIS_FN}: about to delete previously output netCDF file='${TG_NETCDF_FP}'"
-  for CMD in \
-    "rm ${TG_NETCDF_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# MOZART-derived file with source datavar=N2O
-SG_N2O_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_N2O.nc'
-SG_N2O_FN="$(basename ${SG_N2O_URI})"
-export SG_N2O_FP="${SG_DIR}/${SG_N2O_FN}"
-if [[ ! -r "${SG_N2O_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${SG_N2O_FP} ${SG_N2O_URI}" \
-    "ncdump -h ${SG_N2O_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# MOZART-derived file with source datavar=PS
-SG_PS_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
-SG_PS_FN="$(basename ${SG_PS_URI})"
-export SG_PS_FP="${SG_DIR}/${SG_PS_FN}"
-if [[ ! -r "${SG_PS_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${SG_PS_FP} ${SG_PS_URI}" \
-    "ncdump -h ${SG_PS_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# ETOPO1-derived file of surface elevations (only)
-SG_ETOPO1_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
-SG_ETOPO1_FN="$(basename ${SG_ETOPO1_URI})"
-export SG_ETOPO1_FP="${SG_DIR}/${SG_ETOPO1_FN}"
-if [[ ! -r "${SG_ETOPO1_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${SG_ETOPO1_FP} ${SG_ETOPO1_URI}" \
-    "ncdump -h ${SG_ETOPO1_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# MOZART-derived file with all other datavars
-SG_OTHER_URI="${SG_ORIGINAL_URI}"
-SG_OTHER_FN="$(basename ${SG_OTHER_URI})"
-export SG_OTHER_FP="${SG_DIR}/${SG_OTHER_FN}"
-if [[ ! -r "${SG_OTHER_FP}" ]] ; then
-  for CMD in \
-    "wget --no-check-certificate -c -O ${SG_OTHER_FP} ${SG_OTHER_URI}" \
-    "ncdump -h ${SG_OTHER_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-# call NCL script (the real work)---------------------------------------
-
-# ${NCL_EXEC} # bail to NCL and copy script lines
-# TODO: pass commandline args to NCL
-# punt: just use envvars :-(
-for CMD in \
-  "${NCL_EXEC} ${CALL_NCL_FP}" \
-; do
-  echo -e "$ ${CMD}"
-  eval "${CMD}"
-done
-
-# After exiting NCL, show output files ...
-
-if [[ -r "${TG_ASCII_FP}" ||  -r "${TG_NETCDF_FP}" ]] ; then
-  for CMD in \
-    "ls -alht ${WORK_DIR} | head" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-fi
-
-if [[ -r "${TG_ASCII_FP}" ]] ; then
-  for CMD in \
-    "head ${TG_ASCII_FP}" \
-    "tail ${TG_ASCII_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-else # ![[ -r "${TG_ASCII_FP}" ]]
-    echo -e "ERROR: ${THIS_FN}: output ASCII='${TG_ASCII_FP}' not readable"
-fi
-
-if [[ -r "${TG_NETCDF_FP}" ]] ; then
-  for CMD in \
-    "ncdump -h ${TG_NETCDF_FP}" \
-  ; do
-    echo -e "$ ${CMD}"
-    eval "${CMD}"
-  done
-
-# # TODO: plot output
-#   # ... and display output PDF.
-#   if [[ -r "${OUT_PDF_FP}" ]] ; then
-#     for CMD in \
-#       "${PDF_VIEWER} ${OUT_PDF_FP} &" \
-#     ; do
-#       echo -e "$ ${CMD}"
-#       eval "${CMD}"
-#     done
-#   else
-#     echo -e "ERROR: ${THIS_FN}: output PDF='${OUT_PDF_FP}' not readable"
-#   fi
-
-else # ![[ -r "${TG_NETCDF_FP}" ]]
-    echo -e "ERROR: ${THIS_FN}: output netCDF='${TG_NETCDF_FP}' not readable"
-fi

File hsp_to_gh_amsl__MOZART.ncl

View file
+;!/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.
+;----------------------------------------------------------------------
+
+;----------------------------------------------------------------------
+; 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 ; hsp_to_gh_amsl.ncl

File hsp_to_gh_amsl__MOZART.sh

View file
+#!/usr/bin/env bash
+# Requires NCL.
+# Driver for hsp_to_gh_amsl.ncl--configure as needed for your platform.
+# Convert MOZART hybrid sigma-pressure grid (dimensions=(hsp, lat, lon))
+# to a more "Cartesian" grid with z=geometric height above mean sea level
+# (dimensions=(gh_AMSL, lat, lon)). Should work with older NCL versions.
+
+# constants with some simple manipulations------------------------------
+
+THIS_FN="$(basename $0)"
+# will do the real work
+
+# TODO: create from THIS_FN, e.g., use `sed -e s/_/./g`
+export CALL_NCL_FN='hsp_to_gh_amsl.ncl'
+
+# modules on EMVL. TODO: fixme! (these fail from here, are currently just reminders.)
+# You may not need them at all!
+module add netcdf-4.1.2
+#module add ncl_ncarg-6.0.0
+# Temporarily using version on terrae installed by Jerry Herwehe (thanks!)
+#module add nco-4.0.5
+
+# Ensure $NCARG_ROOT is set for NCL--it should be in your dotfiles, but, JIC
+# tlrPanP5
+#NCL_VERSION='ncl_ncarg-6.1.0.Linux_Debian_x86_64_nodap_gcc445'
+#export NCARG_ROOT="${HOME}/bin/${NCL_VERSION}"
+## Only change default if absolutely necessary!
+#NCL_EXEC='ncl'
+# terrae temporary (they gotta upgrade RHEL)
+NCL_VERSION='ncl_ncarg/ncl_ncarg-6.1.0.Linux_RedHat_x86_64_gcc412'
+export NCARG_ROOT="/home/hhg/${NCL_VERSION}"
+NCL_EXEC="${NCARG_ROOT}/bin/ncl"
+
+# workspace
+export WORK_DIR="$(pwd)" # keep it simple for now: same dir as top of repo
+# will do the real work (restricting and plotting)
+export CALL_NCL_FP="${WORK_DIR}/${CALL_NCL_FN}"
+
+# # for plotting
+# PDF_VIEWER='evince'  # whatever works on your platform
+# # temporally disaggregate multiple plots
+# DATE_FORMAT='%Y%m%d_%H%M'
+# export PDF_DIR="${WORK_DIR}"
+
+# source grids
+
+# directories
+SG_DIR="${WORK_DIR}" # keeping it simple
+TG_DIR="${WORK_DIR}"
+
+# files and URIs (from bitbucket repo)
+SG_ORIGINAL_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart.nc'                       # uncropped
+SG_ORIGINAL_FN="$(basename ${SG_ORIGINAL_URI})"
+SG_ORIGINAL_FN_ROOT="${SG_ORIGINAL_FN%.*}"
+SG_ORIGINAL_FN_EXT="${SG_ORIGINAL_FN##*.}"
+NETCDF_EXT="${SG_ORIGINAL_FN_EXT}"                 # stay consistent
+ASCII_EXT='txt'
+
+SG_CROPPED_FN_ROOT="${SG_ORIGINAL_FN_ROOT}_region" # cropped
+SG_CROPPED_FN="${SG_CROPPED_FN_ROOT}.${NETCDF_EXT}"
+
+# MOZART-derived file with source datavar=N2O
+SG_N2O_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_N2O.nc'
+SG_N2O_FN="$(basename ${SG_N2O_URI})"
+export SG_N2O_FP="${SG_DIR}/${SG_N2O_FN}"
+
+# MOZART-derived file with source datavar=PS
+SG_PS_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
+SG_PS_FN="$(basename ${SG_PS_URI})"
+export SG_PS_FP="${SG_DIR}/${SG_PS_FN}"
+
+# ETOPO1-derived file of surface elevations (only)
+SG_ETOPO1_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
+SG_ETOPO1_FN="$(basename ${SG_ETOPO1_URI})"
+export SG_ETOPO1_FP="${SG_DIR}/${SG_ETOPO1_FN}"
+
+# MOZART-derived file with all other datavars
+SG_OTHER_URI="${SG_ORIGINAL_URI}"
+export SG_OTHER_FN="$(basename ${SG_OTHER_URI})" # for doc strings
+export SG_OTHER_FP="${SG_DIR}/${SG_OTHER_FN}"
+
+# target grids
+
+# datavar for elevations of layer interfaces, but horizontally geographical
+export DATAVAR_Z_INT_NAME='z_int_geo'
+# datavar for elevations of layer midpoints, but [ditto]
+export DATAVAR_Z_MID_NAME='z_mid_geo'
+# TODO: (RSN) create fully Cartesian, horizontally Cartesian datavars ...
+
+# target netCDF file ~= ${SG_OTHER_FN} + new datavars above
+TG_NETCDF_FN="${SG_CROPPED_FN_ROOT}_z.${NETCDF_EXT}"
+export TG_NETCDF_FP="${TG_DIR}/${TG_NETCDF_FN}"
+
+# target ASCII file: write ${DATAVAR_Z_INT_NAME} for data interchange with R
+TG_ASCII_FN="${SG_CROPPED_FN_ROOT}_z.${ASCII_EXT}"
+export TG_ASCII_FP="${TG_DIR}/${TG_ASCII_FN}"
+# missing value for string output to be read by R
+export TG_ASCII_NA='NA,NA,NA'
+# not actually usable by @#$%^&! NCL::sprintf
+# export TG_ASCII_LINE_FORMAT='%f,%f,%f' # sprintf string for CSV line
+export TG_ASCII_DATUM_FORMAT='%f' # sprintf string for each datum in line
+
+# payload---------------------------------------------------------------
+
+# setup-----------------------------------------------------------------
+
+mkdir -p ${WORK_DIR}
+
+# get (or remove) data------------------------------------------------
+
+# NCL will not overwrite target ASCII, so we will. (TODO: backup/rename)
+if [[ -r "${TG_ASCII_FP}" ]] ; then
+  echo -e "WARNING: ${THIS_FN}: about to delete previously output ASCII file='${TG_ASCII_FP}'"
+  for CMD in \
+    "rm ${TG_ASCII_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# NCL will not overwrite target netCDF, so we will. (TODO: backup/rename)
+if [[ -r "${TG_NETCDF_FP}" ]] ; then
+  echo -e "WARNING: ${THIS_FN}: about to delete previously output netCDF file='${TG_NETCDF_FP}'"
+  for CMD in \
+    "rm ${TG_NETCDF_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# MOZART-derived file with source datavar=N2O
+SG_N2O_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_N2O.nc'
+SG_N2O_FN="$(basename ${SG_N2O_URI})"
+export SG_N2O_FP="${SG_DIR}/${SG_N2O_FN}"
+if [[ ! -r "${SG_N2O_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${SG_N2O_FP} ${SG_N2O_URI}" \
+    "ncdump -h ${SG_N2O_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# MOZART-derived file with source datavar=PS
+SG_PS_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/2008N2O_restart_region_PS.nc'
+SG_PS_FN="$(basename ${SG_PS_URI})"
+export SG_PS_FP="${SG_DIR}/${SG_PS_FN}"
+if [[ ! -r "${SG_PS_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${SG_PS_FP} ${SG_PS_URI}" \
+    "ncdump -h ${SG_PS_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# ETOPO1-derived file of surface elevations (only)
+SG_ETOPO1_URI='https://bitbucket.org/tlroche/mozart-global-to-aqmeii-na/downloads/ETOPO1_Ice_g_gmt4.grd_region_zS.nc'
+SG_ETOPO1_FN="$(basename ${SG_ETOPO1_URI})"
+export SG_ETOPO1_FP="${SG_DIR}/${SG_ETOPO1_FN}"
+if [[ ! -r "${SG_ETOPO1_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${SG_ETOPO1_FP} ${SG_ETOPO1_URI}" \
+    "ncdump -h ${SG_ETOPO1_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# MOZART-derived file with all other datavars
+SG_OTHER_URI="${SG_ORIGINAL_URI}"
+SG_OTHER_FN="$(basename ${SG_OTHER_URI})"
+export SG_OTHER_FP="${SG_DIR}/${SG_OTHER_FN}"
+if [[ ! -r "${SG_OTHER_FP}" ]] ; then
+  for CMD in \
+    "wget --no-check-certificate -c -O ${SG_OTHER_FP} ${SG_OTHER_URI}" \
+    "ncdump -h ${SG_OTHER_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+# call NCL script (the real work)---------------------------------------
+
+# ${NCL_EXEC} # bail to NCL and copy script lines
+# TODO: pass commandline args to NCL
+# punt: just use envvars :-(
+for CMD in \
+  "${NCL_EXEC} ${CALL_NCL_FP}" \
+; do
+  echo -e "$ ${CMD}"
+  eval "${CMD}"
+done
+
+# After exiting NCL, show output files ...
+
+if [[ -r "${TG_ASCII_FP}" ||  -r "${TG_NETCDF_FP}" ]] ; then
+  for CMD in \
+    "ls -alht ${WORK_DIR} | head" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+fi
+
+if [[ -r "${TG_ASCII_FP}" ]] ; then
+  for CMD in \
+    "head ${TG_ASCII_FP}" \
+    "tail ${TG_ASCII_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+else # ![[ -r "${TG_ASCII_FP}" ]]
+    echo -e "ERROR: ${THIS_FN}: output ASCII='${TG_ASCII_FP}' not readable"
+fi
+
+if [[ -r "${TG_NETCDF_FP}" ]] ; then
+  for CMD in \
+    "ncdump -h ${TG_NETCDF_FP}" \
+  ; do
+    echo -e "$ ${CMD}"
+    eval "${CMD}"
+  done
+
+# # TODO: plot output
+#   # ... and display output PDF.
+#   if [[ -r "${OUT_PDF_FP}" ]] ; then
+#     for CMD in \
+#       "${PDF_VIEWER} ${OUT_PDF_FP} &" \
+#     ; do
+#       echo -e "$ ${CMD}"
+#       eval "${CMD}"
+#     done
+#   else
+#     echo -e "ERROR: ${THIS_FN}: output PDF='${OUT_PDF_FP}' not readable"
+#   fi
+
+else # ![[ -r "${TG_NETCDF_FP}" ]]
+    echo -e "ERROR: ${THIS_FN}: output netCDF='${TG_NETCDF_FP}' not readable"
+fi