"""One chunk of a round 4 historical run. Starts from the L3 gdirs of one v5 experiment and runs them with the climate already in the gdirs: the lmr_seas_fullsplit_oldBE HIST ens183 forcing they were calibrated with. Do not run this directly - see run_slurm.sh and submit.sh. All the settings come from environment variables: R4_GDIR_URL base url of the v5 experiment (.../exps//) R4_PARAMS_FILE the experiment's params.cfg, the same the gdirs were made with R4_EVOLUTION_MODEL empty, or e.g. FluxBased (centerlines) R4_CHUNK_IDX, R4_CHUNK_SIZE R4_RUN_ID, R4_YS R4_DYN_CALIB "1": first calibrate melt_f per glacier with a dynamic run (see below), then do the run with the calibrated melt_f R4_DYN_YS start year of the calibration runs R4_DYN_ERR_FACTOR the Hugonnet error is scaled by this (0.2 as oggm_prepro) R4_DYN_MAXITER max calibration iterations per glacier R4_DYN_STEP_MIN melt_f_max_step_length_minimum (see round4/common.sh) 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 R4_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['R4_GDIR_URL'] PARAMS_FILE = os.environ['R4_PARAMS_FILE'] EVOLUTION_MODEL = os.environ.get('R4_EVOLUTION_MODEL', '') CHUNK_IDX = int(os.environ['R4_CHUNK_IDX']) CHUNK_SIZE = int(os.environ['R4_CHUNK_SIZE']) RUN_ID = os.environ['R4_RUN_ID'] YS = int(os.environ['R4_YS']) DYN_CALIB = os.environ.get('R4_DYN_CALIB', '') == '1' DYN_YS = int(os.environ.get('R4_DYN_YS', 1000)) DYN_ERR_FACTOR = float(os.environ.get('R4_DYN_ERR_FACTOR', 0.2)) DYN_MAXITER = int(os.environ.get('R4_DYN_MAXITER', 20)) DYN_STEP_MIN = float(os.environ.get('R4_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. # Not stored in the gdirs: oggm_prepro sets it at run time for centerlines if EVOLUTION_MODEL: cfg.PARAMS['evolution_model'] = EVOLUTION_MODEL 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'] 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'] # 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)