adding option for constant lapserate

pygem/pygem_input.py

@@ -236,6 +236,7 @@
 
 #%% ===== MODEL PARAMETERS =====
 use_calibrated_modelparams = True   # False: use input values, True: use calibrated model parameters
+use_constant_lapserate = False      # False: use spatially and temporally varying lapse rate, True: use constant value specified below
 if not use_calibrated_modelparams:
     print('\nWARNING: using non-calibrated model parameters\n')
     sim_iters = 1
@@ -246,8 +247,7 @@
 ddfsnow_iceratio = 0.7              # Ratio degree-day factor snow snow to ice
 ddfice = ddfsnow / ddfsnow_iceratio # degree-day factor of ice [m w.e. d-1 degC-1]
 precgrad = 0.0001                   # precipitation gradient on glacier [m-1]
-lrgcm = -0.0065                     # lapse rate from gcm to glacier [K m-1]
-lrglac = -0.0065                    # lapse rate on glacier for bins [K m-1]
+lapserate = -0.0065                 # temperature lapse rate for both gcm to glacier and on glacier between elevation bins [K m-1]
 tsnow_threshold = 1                 # temperature threshold for snow [deg C] (HH2015 used 1.5 degC +/- 1 degC)
 calving_k = 0.7                     # frontal ablation rate [yr-1]
 
@@ -318,7 +318,8 @@
     assert os.path.exists(era5_fp + era5_temp_fn), 'ERA5 temperature filepath does not exist'
     assert os.path.exists(era5_fp + era5_prec_fn), 'ERA5 precipitation filepath does not exist'
     assert os.path.exists(era5_fp + era5_elev_fn), 'ERA5 elevation data does not exist'
-    assert os.path.exists(era5_fp + era5_lr_fn), 'ERA5 lapse rate data does not exist'
+    if not use_constant_lapserate:
+        assert os.path.exists(era5_fp + era5_lr_fn), 'ERA5 lapse rate data does not exist'
     if option_ablation == 2:
         assert os.path.exists(era5_fp + era5_tempstd_fn), 'ERA5 temperature std filepath does not exist'
 

run_calibration_frontalablation.py

@@ -122,7 +122,10 @@ def reg_calving_flux(main_glac_rgi, calving_k, fa_glac_data_reg=None,
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
     # Lapse rate [degC m-1]
-    gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+    if pygem_prms.use_constant_lapserate:
+        gcm_lr = np.zeros(gcm_temp.shape) + pygem_prms.lapserate
+    else:
+        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
 
     # ===== CALIBRATE ALL THE GLACIERS AT ONCE =====
     output_cns = ['RGIId', 'calving_k', 'calving_thick', 'calving_flux_Gta_inv', 'calving_flux_Gta', 'no_errors', 'oggm_dynamics']

run_calibration_icethickness_consensus.py

@@ -264,7 +264,10 @@ def reg_vol_comparison(gdirs, mbmods, a_multiplier=1, fs=0, debug=False):
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi_subset)
     # Lapse rate [degC m-1]
-    gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi_subset, dates_table)
+    if pygem_prms.use_constant_lapserate:
+        gcm_lr = np.zeros(gcm_temp.shape) + pygem_prms.lapserate
+    else:
+        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
 
     # ===== RUN MASS BALANCE =====
     # Number of years (for OGGM's run_until_and_store)

run_calibration_woggm.py

@@ -346,7 +346,11 @@ def main(list_packed_vars):
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
     # Lapse rate [degC m-1]
-    gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+    if pygem_prms.use_constant_lapserate:
+        gcm_lr = np.zeros(gcm_temp.shape) + pygem_prms.lapserate
+    else:
+        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+
 
     # ===== LOOP THROUGH GLACIERS TO RUN CALIBRATION =====
     for glac in range(main_glac_rgi.shape[0]):

run_preprocessing_wgms_mbdata.py

@@ -228,8 +228,11 @@ def getparser():
                                                                   dates_table)
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
-    # Lapse rate
-    gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+    # Lapse rate [degC m-1]
+    if pygem_prms.use_constant_lapserate:
+        gcm_lr = np.zeros(gcm_temp.shape) + pygem_prms.lapserate
+    else:
+        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
 
     # ===== PROCESS THE OBSERVATIONS ======
     prec_cn = pygem_prms.ref_gcm_name + '_prec'

run_simulation_woggm.py

@@ -1006,15 +1006,18 @@ def main(list_packed_vars):
                                                                   dates_table)
     # Elevation [m asl]
     gcm_elev = gcm.importGCMfxnearestneighbor_xarray(gcm.elev_fn, gcm.elev_vn, main_glac_rgi)
-    # Lapse rate
-    if gcm_name in ['ERA-Interim', 'ERA5']:
-        gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+    # Lapse rate [degC m-1]
+    if pygem_prms.use_constant_lapserate:
+        gcm_lr = np.zeros(gcm_temp.shape) + pygem_prms.lapserate
     else:
-        # Compute lapse rates based on reference climate data
-        ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi,
-                                                                        dates_table_ref)
-        # Monthly average from reference climate data
-        gcm_lr = gcmbiasadj.monthly_avg_array_rolled(ref_lr, dates_table_ref, dates_table)
+        if gcm_name in ['ERA-Interim', 'ERA5']:
+            gcm_lr, gcm_dates = gcm.importGCMvarnearestneighbor_xarray(gcm.lr_fn, gcm.lr_vn, main_glac_rgi, dates_table)
+        else:
+            # Compute lapse rates based on reference climate data
+            ref_lr, ref_dates = ref_gcm.importGCMvarnearestneighbor_xarray(ref_gcm.lr_fn, ref_gcm.lr_vn, main_glac_rgi,
+                                                                            dates_table_ref)
+            # Monthly average from reference climate data
+            gcm_lr = gcmbiasadj.monthly_avg_array_rolled(ref_lr, dates_table_ref, dates_table)
 
 
     # ===== BIAS CORRECTIONS =====