"""One chunk of a round 4a run, from the gdirs of one v6 experiment. Copy of round4/run_hist.py (V6_ instead of R4_ variables, no centerlines, stop reason in the table). Starts from the L3 gdirs of one v6 experiment and runs them with the climate already in the gdirs (the forcing they were calibrated with). Do not run this directly - see run_dyn_slurm.sh and submit.sh. All the settings come from environment variables: V6_GDIR_URL base url of the v6 experiment (.../exps//) V6_PARAMS_FILE the experiment's params.cfg, the same the gdirs were made with V6_CHUNK_IDX, V6_CHUNK_SIZE V6_RUN_ID, V6_YS V6_DYN_CALIB "1": first calibrate melt_f per glacier with a dynamic run (see below), then do the run with the calibrated melt_f V6_DYN_YS start year of the calibration runs V6_DYN_ERR_FACTOR the Hugonnet error is scaled by this (0.2 as oggm_prepro) V6_DYN_MAXITER max calibration iterations per glacier V6_DYN_STEP_MIN melt_f_max_step_length_minimum (see common.sh) V6_TEST_IDS optional: only these glaciers (smoke test) OGGM_WORKDIR, OGGM_OUTDIR The dynamic calibration is OGGM's run_dynamic_melt_f_calibration with the plain dynamic_melt_f_run: no dynamic spinup, i.e. no temp_bias loop matching the RGI area. Each iteration runs the glacier from the RGI-date flowlines in V6_DYN_YS to the end of the climate, and compares its 2000-2020 mass balance with the Hugonnet value in the gdir (observations.yml). Only melt_f changes: prcp_fac and temp_bias keep their static calibration values. """ import os import logging import geopandas as gpd import pandas as pd import oggm.cfg as cfg from oggm import utils, workflow, tasks GDIR_URL = os.environ['V6_GDIR_URL'] PARAMS_FILE = os.environ['V6_PARAMS_FILE'] CHUNK_IDX = int(os.environ['V6_CHUNK_IDX']) CHUNK_SIZE = int(os.environ['V6_CHUNK_SIZE']) RUN_ID = os.environ['V6_RUN_ID'] YS = int(os.environ['V6_YS']) DYN_CALIB = os.environ.get('V6_DYN_CALIB', '') == '1' DYN_YS = int(os.environ.get('V6_DYN_YS', 1000)) DYN_ERR_FACTOR = float(os.environ.get('V6_DYN_ERR_FACTOR', 0.2)) DYN_MAXITER = int(os.environ.get('V6_DYN_MAXITER', 20)) DYN_STEP_MIN = float(os.environ.get('V6_DYN_STEP_MIN', 0.1)) RGI_REG = '11' # Same params as when the gdirs were made. The calibrated values (melt_f, # prcp_fac, temp_bias, glen_a...) live in each gdir's settings.yml and take # precedence over these anyway. cfg.initialize(file=PARAMS_FILE, logging_level='WORKFLOW') cfg.PARAMS['border'] = 160 WORKING_DIR = os.environ.get('OGGM_WORKDIR', '') if not WORKING_DIR: raise RuntimeError('Need a working dir') utils.mkdir(WORKING_DIR) cfg.PATHS['working_dir'] = WORKING_DIR OUTPUT_DIR = os.environ.get('OGGM_OUTDIR', '') if not OUTPUT_DIR: raise RuntimeError('Need an output dir') utils.mkdir(OUTPUT_DIR) cfg.PARAMS['use_multiprocessing'] = True cfg.PARAMS['continue_on_error'] = True # round3's variables, plus mass: the dynamic calibration needs it cfg.PARAMS['store_diagnostic_variables'] = ['volume', 'volume_bsl', 'area', 'area_min_h', 'length', 'mass'] cfg.PARAMS['min_ice_thick_for_area'] = 2. cfg.PARAMS['min_ice_thick_for_length'] = 2. workflow.init_mp_pool(True) log = logging.getLogger(__name__) # Our chunk of glaciers: the same blocks of the RGI id space as the chunked # prepro (see workflow.get_rgi_chunk) rgidf = gpd.read_file(utils.get_rgi_region_file(RGI_REG, version='62')) rgidf = workflow.get_rgi_chunk(rgidf, CHUNK_IDX, chunk_size=CHUNK_SIZE) # The Alps only: 11-02 (Pyrenees) is outside the LMR grid, and its gdirs fail # the mass balance calibration anyway rgidf = rgidf.loc[rgidf.O2Region == '1'] TEST_IDS = os.environ.get('V6_TEST_IDS', '').split() if TEST_IDS: rgidf = rgidf.loc[rgidf.RGIId.isin(TEST_IDS)] if len(rgidf) == 0: # The ids have gaps, so an empty chunk is a normal thing log.workflow(f'Chunk {CHUNK_IDX} is empty, nothing to do.') raise SystemExit(0) log.workflow(f'Starting {RUN_ID} for RGI reg {RGI_REG}, chunk {CHUNK_IDX} ' f'(size {CHUNK_SIZE}): {len(rgidf)} glaciers, ' f'evolution model {cfg.PARAMS["evolution_model"]}') gdirs = workflow.init_glacier_directories(rgidf, from_prepro_level=3, prepro_base_url=GDIR_URL, prepro_rgi_version='62', prepro_border=160) if DYN_CALIB: from oggm.core.dynamic_spinup import (run_dynamic_melt_f_calibration, dynamic_melt_f_run, dynamic_melt_f_run_fallback) log.workflow(f'Dynamic melt_f calibration, chunk {CHUNK_IDX}: runs from ' f'{DYN_YS}, Hugonnet error x {DYN_ERR_FACTOR}, ' f'maxiter {DYN_MAXITER}, min step {DYN_STEP_MIN}') workflow.execute_entity_task(run_dynamic_melt_f_calibration, gdirs, ref_mb_err_scaling_factor=DYN_ERR_FACTOR, maxiter=DYN_MAXITER, melt_f_max_step_length_minimum=DYN_STEP_MIN, ys=DYN_YS, run_function=dynamic_melt_f_run, fallback_function=dynamic_melt_f_run_fallback, ignore_errors=True, output_filesuffix='_dyn_calib') # What the calibration did, per glacier keys = ['used_spinup_option', 'melt_f_before_dynamic_calibration', 'melt_f_dynamic_calibration', 'melt_f', 'dmdtda_mismatch_dynamic_calibration_reference', 'dmdtda_dynamic_calibration_given_error', 'dmdtda_mismatch_with_initial_melt_f', 'dmdtda_mismatch_dynamic_calibration', 'run_dynamic_melt_f_calibration_iterations', 'run_dynamic_melt_f_calibration_stop_reason'] # Read from diagnostics.json: the calibration ran in the worker processes, # and the settings object of the gdirs in this process is a stale copy rows = {} for gdir in gdirs: diag = gdir.get_diagnostics() r = {k: diag.get(k) for k in keys} r['melt_f'] = diag.get('melt_f_dynamic_calibration') r['error'] = gdir.get_error_log() rows[gdir.rgi_id] = r pd.DataFrame(rows).T.to_csv(os.path.join( OUTPUT_DIR, f'dyn_calib_reg{RGI_REG}_chunk{CHUNK_IDX:03d}.csv')) # The run, with the calibrated melt_f if any (it is in the gdir settings) workflow.execute_entity_task(tasks.run_from_climate_data, gdirs, climate_filename='climate_historical', output_filesuffix=RUN_ID, return_value=False, ys=YS) log.workflow(f'Aggregating {RUN_ID}, chunk {CHUNK_IDX}') out_nc_path = os.path.join(OUTPUT_DIR, f'{RUN_ID}_reg{RGI_REG}_chunk{CHUNK_IDX:03d}.nc') utils.compile_run_output(gdirs, input_filesuffix=RUN_ID, path=out_nc_path)