Allow for frontal ablation in dynamical spinup (#175)

* Allow for calving in dynamical spinup
* Pull cfl params from pygem config
* Allow for inversion in mcmc elev_change_1d calib
* Optimize spinup run 0-year period with inversion dynamical run
* Store model geom
* Bug fix in indexing of time steps in mass balance loop introduced when updating to allow for daily calibration

Author: btobers

poetry.lock

@@ -156,13 +156,19 @@ def run(simpath, debug=False):
             gdir,
             fl_diag=pygem_fl_diag,
             concat_input_filesuffix='_spinup_historical',  # concatenate with the historical spinup
-            output_filesuffix=f'_pygem_{f_suffix}',  # filesuffix added to the output filename gridded_simulation.nc, if empty input_filesuffix is used
+            output_filesuffix=f'_{f_suffix}',  # filesuffix added to the output filename gridded_simulation.nc, if empty input_filesuffix is used
         )[0]
-        print(
-            '2D simulated ice thickness created: ',
-            gdir.get_filepath('gridded_simulation', filesuffix=f'_pygem_{f_suffix}'),
-        )
+
+        # copy glacier_ext from gridded_data
+        with xr.open_dataset(gdir.get_filepath('gridded_data')) as ds0:
+            ds['glacier_ext'] = ds0['glacier_ext']
+            ds.to_netcdf(gdir.get_filepath('gridded_simulation', filesuffix=f'_{f_suffix}'))
+
         if debug:
+            print(
+                '2D simulated ice thickness created: ',
+                gdir.get_filepath('gridded_simulation', filesuffix=f'_{f_suffix}'),
+            )
             plot_distributed_thickness(ds)
 
     return

pygem/bin/postproc/postproc_distribute_ice.py

@@ -354,13 +354,10 @@ def calc_thick_change_1d(gdir, fls, mbmod, ds):
     years_subannual = np.array([d.year for d in gdir.dates_table['date']])
     yrs = np.unique(years_subannual)
     # grab components of interest
-    bin_thick_annual = ds[0].thickness_m.values.T  # glacier thickness [m ice], (nbins, nyears)
-    bin_massbalclim_ice_annual = ds[0].climatic_mb_myr.values.T[
+    bin_massbalclim_ice_annual = ds[0].climatic_mb.values.T[
         :, 1:
     ]  # annual climatic mass balance [m ice] (nbins, nyears - 1)
-    bin_delta_thick_annual = ds[0].dhdt_myr.values.T[
-        :, 1:
-    ]  # annual change in ice thickness [m ice] (nbins, nyears - 1)
+    bin_delta_thick_annual = ds[0].dhdt.values.T[:, 1:]  # annual change in ice thickness [m ice] (nbins, nyears - 1)
     bin_flux_divergence_annual = (
         bin_massbalclim_ice_annual - bin_delta_thick_annual
     )  # annual flux divergence [m ice] (nbins, nyears - 1)
@@ -384,8 +381,9 @@ def calc_thick_change_1d(gdir, fls, mbmod, ds):
     bin_delta_thick_subannual = bin_massbalclim_ice - bin_flux_divergence_subannual
 
     # --- Step 4: calculate subannual thickness = running thickness change + initial thickness---
+    bin_thick_initial = ds[0].thickness_m.isel(time=0).values  # initial glacier thickness [m ice], (nbins)
     running_bin_delta_thick_subannual = np.nancumsum(bin_delta_thick_subannual, axis=-1)
-    bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_annual[:, 0][:, np.newaxis]
+    bin_thick_subannual = running_bin_delta_thick_subannual + bin_thick_initial[:, np.newaxis]
 
     # --- Step 5: rebin ---
     # get surface height at the specified reference year

pygem/bin/run/run_calibration.py

@@ -401,13 +401,11 @@ def main():
         nargs='+',
         help='Randoph Glacier Inventory glacier number (can take multiple)',
     )
-    (
-        parser.add_argument(
-            '-rgi_glac_number_fn',
-            type=str,
-            default=None,
-            help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
-        ),
+    parser.add_argument(
+        '-rgi_glac_number_fn',
+        type=str,
+        default=None,
+        help='Filepath containing list of rgi_glac_number, helpful for running batches on spc',
     )
     parser.add_argument(
         '-calibrate_regional_glen_a',
@@ -473,10 +471,10 @@ def main():
             # format appropriately
             glac_no = [float(g) for g in glac_no]
             batches = [f'{g:.5f}' if g >= 10 else f'0{g:.5f}' for g in glac_no]
-            regional_inv = False  # flag to indicate per-glacier inversion
         elif args.rgi_glac_number_fn is not None:
             with open(args.rgi_glac_number_fn, 'r') as f:
                 batches = json.load(f)
+        regional_inv = False  # flag to indicate per-glacier inversion
 
     # set up partial function with common arguments
     run_partial = partial(