MOM support

.gitignore

@@ -105,3 +105,6 @@ venv.bak/
 */logs/*
 */dask-worker-space/*
 */notes/*
+
+# netCDF output
+**/*.nc

Fe_sediment_flux/sedfrac_compute.ipynb

@@ -652,9 +652,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python [conda env:dev]",
+   "display_name": "Python [conda env:miniconda3-python-tutorial]",
    "language": "python",
-   "name": "conda-env-dev-py"
+   "name": "conda-env-miniconda3-python-tutorial-py"
   },
   "language_info": {
    "codemirror_mode": {
@@ -670,5 +670,5 @@
   }
  },
  "nbformat": 4,
- "nbformat_minor": 2
+ "nbformat_minor": 4
 }

abio/gen_files.sh

@@ -0,0 +1,32 @@
+#!/bin/bash -e
+
+DATE=`date "+%y%m%d"`
+ORIG_FILE=/glade/campaign/cesm/cesmdata/inputdata/lnd/clm2/isotopes/atm_delta_C14_CMIP6_3x1_global_1850-2015_yearly_v2.0_c190528.nc
+OUTPUT_DIR=/glade/work/mlevy/cesm_inputdata
+
+rm -f ../atm_delta_C14_CMIP6*.nc
+
+for sector in `seq 1 3`
+do
+  OUT_FILE=${OUTPUT_DIR}/atm_delta_C14_CMIP6_sector${sector}_global_1850-2015_yearly_v2.0_c${DATE}.nc
+  # (1) pull out sector number ${sector}
+  ncks -F -d sector,${sector} ${ORIG_FILE} ${OUT_FILE}
+  # (2) average over single sector (i.e. remove sector dim from data)
+  ncwa -O -a sector ${OUT_FILE} ${OUT_FILE}
+  # (3) make time an unlimited dimesnion
+  ncks -O --mk_rec_dmn time -x -v sector ${OUT_FILE} ${OUT_FILE}
+  # (4) pull out first and last time levels to get Jan 1, 1850 and Jan 1, 2016 dates
+  ncks -F -d time,1 ${OUT_FILE} .tmp_pre.nc
+  ncap2 -O -F -s 'time(1)=0.0' .tmp_pre.nc .tmp_pre.nc
+  ncks -F -d time,166 ${OUT_FILE} .tmp_post.nc
+  ncap2 -O -F -s 'time(1)=166.0' .tmp_post.nc .tmp_post.nc
+  # (5) concatenate all the files
+  ncrcat -O .tmp_pre.nc ${OUT_FILE} .tmp_post.nc ${OUT_FILE}
+  # (6) update time bounds
+  ncap2 -O -F -s 'bounds_time(1,2)=0.5;bounds_time(2,1)=0.5' ${OUT_FILE} ${OUT_FILE}
+  ncap2 -O -F -s 'bounds_time(167,2)=165.5;bounds_time(168,1)=165.5' ${OUT_FILE} ${OUT_FILE}
+  rm -f .tmp_pre.nc .tmp_post.nc
+
+  echo "Done processing sector ${sector}:"
+  ncdump -v bounds_sector ${OUT_FILE} | grep "bounds_sector ="
+done

initial_conditions/MOM_grid.yaml

@@ -0,0 +1,15 @@
+lateral_dims: [ 458,  540]
+vertical_dims: 34
+type: curvilinear
+open_vertical_grid:
+    vert_grid_file: grid_data/tx0.66v1_vert_grid
+    depth_units: m
+    depth_coord_name: z_t
+gen_grid_file:
+    method: gen_POP_grid_file
+    kwargs:
+        horiz_grid_fname: /glade/p/cesmdata/cseg/inputdata/ocn/mom/tx0.66v1/grid/horiz_grid_20190315.ieeer8
+        topography_fname: /glade/p/cesmdata/cseg/inputdata/ocn/mom/tx0.66v1/grid/topography_20190315.ieeei4
+        region_mask_fname: None
+        type: tripole
+        lateral_dims: [ 458,  540]

ndep/MOM_grid.yaml

@@ -0,0 +1,15 @@
+lateral_dims: [ 458,  540]
+vertical_dims: 34
+type: curvilinear
+open_vertical_grid:
+    vert_grid_file: grid_data/tx0.66v1_vert_grid
+    depth_units: m
+    depth_coord_name: z_t
+gen_grid_file:
+    method: gen_POP_grid_file
+    kwargs:
+        horiz_grid_fname: /glade/p/cesmdata/cseg/inputdata/ocn/mom/tx0.66v1/grid/horiz_grid_20190315.ieeer8
+        topography_fname: /glade/p/cesmdata/cseg/inputdata/ocn/mom/tx0.66v1/grid/topography_20190315.ieeei4
+        region_mask_fname: None
+        type: tripole
+        lateral_dims: [ 458,  540]

ndep/config.py

@@ -0,0 +1,13 @@
+import os
+
+USER = os.environ['USER']
+
+DATAROOT = f'/glade/work/{USER}'
+
+inputdata = '/glade/p/cesmdata/cseg/inputdata'
+
+dirwork = f'{DATAROOT}/cesm_inputdata/work'
+os.makedirs(dirwork, exist_ok=True)   
+
+dirout = f'{DATAROOT}/cesm_inputdata'
+os.makedirs(dirout, exist_ok=True)   

ndep/gen_remap_weights_1d.py

@@ -0,0 +1,42 @@
+import cf_units
+import numpy as np
+from scipy.sparse import csr_matrix
+
+def interval_overlap(interval_1, interval_2):
+    '''
+    compute how much overlap there is between two intervals
+
+    intervals can be tuples, lists, numpy arrays
+    intervals can be specified in an increasing or decreasing order
+    '''
+    overlap = min([max(interval_1), max(interval_2)]) - max([min(interval_1), min(interval_2)])
+    return max([0.0, overlap])
+
+
+def gen_remap_weights_1d(src_bnds, dst_bnds, src_units='m', dst_units='m'):
+    '''
+    generate weights to remap from 1d src axis to 1d dst axis
+
+    axes are represented as interval bounds, provided in a 2D xarray DataArray
+
+    weights are returned as a scipy.sparse matrix
+    '''
+
+    src_units_obj = cf_units.Unit(src_units)
+    dst_units_obj = cf_units.Unit(dst_units)
+    src_len = src_bnds.shape[0]
+    dst_len = dst_bnds.shape[0]
+    weights_dense = np.zeros((dst_len, src_len))
+    for src_ind in range(src_len):
+        src_values_dst_units = src_units_obj.convert(src_bnds[src_ind, :].values, dst_units_obj)
+        for dst_ind in range(dst_len):
+            weights_dense[dst_ind, src_ind] = \
+                interval_overlap(src_values_dst_units, dst_bnds[dst_ind, :].values)
+
+    # normalize rows to have sum 1.0
+    for dst_ind in range(dst_len):
+        row_sum = weights_dense[dst_ind, :].sum()
+        if row_sum > 0.0:
+            weights_dense[dst_ind, :] /= row_sum
+
+    return csr_matrix(weights_dense)