From 8cbcd22636d4362a01838655b1794901550b05dc Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Wed, 6 May 2026 18:15:53 -0700 Subject: [PATCH 01/15] Add PROTECT emulandice2 experiment for reference. Update emulandice2 pipelines to have build files for Mar2024 PROTECT deliverable --- .../emulandice2.protect.ssp585/config.yml | 26 ++++ modules/emulandice2/AIS/pipeline.2300.yml | 66 +++++++++ modules/emulandice2/AIS/pipeline.yml | 4 +- modules/emulandice2/GrIS/pipeline.2300.yml | 66 +++++++++ modules/emulandice2/GrIS/pipeline.yml | 4 +- .../emulandice2/glaciers/pipeline.2300.yml | 137 ++++++++++++++++++ modules/emulandice2/glaciers/pipeline.yml | 2 +- 7 files changed, 300 insertions(+), 5 deletions(-) create mode 100644 experiments/emulandice2.protect.ssp585/config.yml create mode 100644 modules/emulandice2/AIS/pipeline.2300.yml create mode 100644 modules/emulandice2/GrIS/pipeline.2300.yml create mode 100644 modules/emulandice2/glaciers/pipeline.2300.yml diff --git a/experiments/emulandice2.protect.ssp585/config.yml b/experiments/emulandice2.protect.ssp585/config.yml new file mode 100644 index 00000000..a8baf4a8 --- /dev/null +++ b/experiments/emulandice2.protect.ssp585/config.yml @@ -0,0 +1,26 @@ +global-options: + nsamps: 2237 + scenario: ssp585 + pyear_start: 2020 + pyear_end: 2300 #2150 + pyear_step: 10 + baseyear: 2005 + pipeline_file: "pipeline.2300.yml" # "pipeline.yml" + +climate_step: + temperature: + module_set: "fair" + module: "temperature" + pipeline_file: "pipeline.global.yml" + generates_climate_output: true + seed: 2024 +sealevel_step: + emuAIS: + module_set: "emulandice2" + module: "AIS" + emuGLA: + module_set: "emulandice2" + module: "glaciers" + emuGIS: + module_set: "emulandice2" + module: "GrIS" diff --git a/modules/emulandice2/AIS/pipeline.2300.yml b/modules/emulandice2/AIS/pipeline.2300.yml new file mode 100644 index 00000000..49a4cabd --- /dev/null +++ b/modules/emulandice2/AIS/pipeline.2300.yml @@ -0,0 +1,66 @@ + +project: + task1: + python_dependencies: "xarray dask numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - AIS + - --region + - ALL + - --emu_file + - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_AIS_project_2300_240214_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_ALL_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_ALL_globalsl.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/AIS/pipeline.yml b/modules/emulandice2/AIS/pipeline.yml index 50654cb7..c3f6ec22 100644 --- a/modules/emulandice2/AIS/pipeline.yml +++ b/modules/emulandice2/AIS/pipeline.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + - emu_file/AIS_ALL_Kori_PISM_CISM_ElmerIce_pow_exp_10_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -27,7 +27,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_data.tgz" + - "emulandice2_AIS_project_2150_240214_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.2300.yml b/modules/emulandice2/GrIS/pipeline.2300.yml new file mode 100644 index 00000000..39d102a6 --- /dev/null +++ b/modules/emulandice2/GrIS/pipeline.2300.yml @@ -0,0 +1,66 @@ + +project: + task1: + python_dependencies: "xarray dask numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - GIS + - --region + - ALL + - --emu_file + - emu_file/GIS_ALL_CISM_pow_exp_01_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_GrIS_project_2300_240214_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_ALL_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_ALL_globalsl.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/GrIS/pipeline.yml b/modules/emulandice2/GrIS/pipeline.yml index e5cf200a..aa9f0815 100644 --- a/modules/emulandice2/GrIS/pipeline.yml +++ b/modules/emulandice2/GrIS/pipeline.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/GIS_ALL_CISM_pow_exp_01_EMULATOR.RData + - emu_file/GIS_ALL_CISM_IMAUICE_ElmerIce_GISM_pow_exp_01_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -27,7 +27,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_data.tgz" + - "emulandice2_GrIS_project_2100_240214_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/glaciers/pipeline.2300.yml b/modules/emulandice2/glaciers/pipeline.2300.yml new file mode 100644 index 00000000..d479b6b6 --- /dev/null +++ b/modules/emulandice2/glaciers/pipeline.2300.yml @@ -0,0 +1,137 @@ + +project: + task1: + python_dependencies: "xarray dask numpy scipy" + executable: "python3" + script_noupload: "emulandice_project.py" + arguments: + - --ice_source + - GLA + - --region + - RGI01 + - RGI02 + - RGI03 + - RGI04 + - RGI05 + - RGI06 + - RGI07 + - RGI08 + - RGI09 + - RGI10 + - RGI11 + - RGI12 + - RGI13 + - RGI14 + - RGI15 + - RGI16 + - RGI17 + - RGI18 + - RGI19 + - --emu_file + - emu_file/GLA_RGI01_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI02_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI03_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI04_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI05_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI06_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI07_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + - emu_file/GLA_RGI08_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI09_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI10_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI11_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI12_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI13_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI14_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI15_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI16_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI17_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI18_GloGEM_OGGM_pow_exp_20_EMULATOR.RData + - emu_file/GLA_RGI19_GloGEM_OGGM_pow_exp_01_EMULATOR.RData + options: + - "pipeline_id" + - "scenario" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_input_data: + - "%MODULE_PATH%/../shared/emulandice_environment.sh" + - "%MODULE_PATH%/../shared/emulandice_steer.sh" + - "%MODULE_PATH%/../shared/emulandice2/main.R" + - "%MODULE_PATH%/../shared/emulandice_project.py" + upload_and_extract_input_data: + - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" + - "emulandice2_glaciers_project_2300_240223_data.tgz" + global_total_files: + - "%PIPELINE_ID%_ALL_globalsl.nc" + download_output_data: + - "%PIPELINE_ID%_ALL_globalsl.nc" + +postprocess: + task1: + executable: "python3" + script_noupload: "emulandice_postprocess.py" + arguments: + - --ncfiles + - "%PIPELINE_ID%_RGI01_globalsl.nc" + - "%PIPELINE_ID%_RGI02_globalsl.nc" + - "%PIPELINE_ID%_RGI03_globalsl.nc" + - "%PIPELINE_ID%_RGI04_globalsl.nc" + - "%PIPELINE_ID%_RGI05_globalsl.nc" + - "%PIPELINE_ID%_RGI06_globalsl.nc" + - "%PIPELINE_ID%_RGI07_globalsl.nc" + - "%PIPELINE_ID%_RGI08_globalsl.nc" + - "%PIPELINE_ID%_RGI09_globalsl.nc" + - "%PIPELINE_ID%_RGI10_globalsl.nc" + - "%PIPELINE_ID%_RGI11_globalsl.nc" + - "%PIPELINE_ID%_RGI12_globalsl.nc" + - "%PIPELINE_ID%_RGI13_globalsl.nc" + - "%PIPELINE_ID%_RGI14_globalsl.nc" + - "%PIPELINE_ID%_RGI15_globalsl.nc" + - "%PIPELINE_ID%_RGI16_globalsl.nc" + - "%PIPELINE_ID%_RGI17_globalsl.nc" + - "%PIPELINE_ID%_RGI18_globalsl.nc" + - "%PIPELINE_ID%_RGI19_globalsl.nc" + options: + - "locationfile" + - "chunksize" + - "pipeline_id" + - "scenario" + - "baseyear" + - "grdfingerprintfile" + upload_input_data: + - '%MODULE_PATH%/../shared/read_locationfile.py' + - '%MODULE_PATH%/../shared/AssignFP.py' + - '%MODULE_PATH%/../shared/ReadFingerprint.py' + - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' + - '%MODULE_PATH%/../shared/emulandice_postprocess.py' + upload_and_extract_input_data: + - 'grd_fingerprints_data.tgz' + copy_shared_data: + - '$SHARED/location.lst' + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_RGI01_globalsl.nc" + - "%PIPELINE_ID%_RGI02_globalsl.nc" + - "%PIPELINE_ID%_RGI03_globalsl.nc" + - "%PIPELINE_ID%_RGI04_globalsl.nc" + - "%PIPELINE_ID%_RGI05_globalsl.nc" + - "%PIPELINE_ID%_RGI06_globalsl.nc" + - "%PIPELINE_ID%_RGI07_globalsl.nc" + - "%PIPELINE_ID%_RGI08_globalsl.nc" + - "%PIPELINE_ID%_RGI09_globalsl.nc" + - "%PIPELINE_ID%_RGI10_globalsl.nc" + - "%PIPELINE_ID%_RGI11_globalsl.nc" + - "%PIPELINE_ID%_RGI12_globalsl.nc" + - "%PIPELINE_ID%_RGI13_globalsl.nc" + - "%PIPELINE_ID%_RGI14_globalsl.nc" + - "%PIPELINE_ID%_RGI15_globalsl.nc" + - "%PIPELINE_ID%_RGI16_globalsl.nc" + - "%PIPELINE_ID%_RGI17_globalsl.nc" + - "%PIPELINE_ID%_RGI18_globalsl.nc" + - "%PIPELINE_ID%_RGI19_globalsl.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/glaciers/pipeline.yml b/modules/emulandice2/glaciers/pipeline.yml index 4be08ffc..b8f21b30 100644 --- a/modules/emulandice2/glaciers/pipeline.yml +++ b/modules/emulandice2/glaciers/pipeline.yml @@ -63,7 +63,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_glaciers_project_data.tgz" + - "emulandice2_glaciers_project_2100_240223_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: From c92a49ab9fc477d70ad24846accc0682511aea96 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Wed, 6 May 2026 18:33:29 -0700 Subject: [PATCH 02/15] fix emulandice2.protect experiment config --- experiments/emulandice2.protect.ssp585/config.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/experiments/emulandice2.protect.ssp585/config.yml b/experiments/emulandice2.protect.ssp585/config.yml index a8baf4a8..c1cb212f 100644 --- a/experiments/emulandice2.protect.ssp585/config.yml +++ b/experiments/emulandice2.protect.ssp585/config.yml @@ -13,7 +13,6 @@ climate_step: module: "temperature" pipeline_file: "pipeline.global.yml" generates_climate_output: true - seed: 2024 sealevel_step: emuAIS: module_set: "emulandice2" From 54acc90ec56d3ee14ae3c26a519059444b1447dc Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 11 May 2026 10:32:34 -0700 Subject: [PATCH 03/15] allow h5 _climate file --- FACTS.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/FACTS.py b/FACTS.py index d46c9f98..18e7cb64 100644 --- a/FACTS.py +++ b/FACTS.py @@ -395,7 +395,7 @@ def IdentifyClimateOutputFiles(pcfg,pipe_name): tcfg = pcfg[this_stage][this_task] if "climate_output_data" in tcfg.keys(): for this_file in tcfg['climate_output_data']: - if this_file.__contains__('climate.nc'): + if (this_file.__contains__('climate.nc') or this_file.__contains__('climate.h5')): pd['climate'] = '$SHARED/climate/' + mvar_replace_dict(mvar_dict, this_file) elif this_file.__contains__('gsat.nc'): pd['gsat'] = '$SHARED/climate/' + mvar_replace_dict(mvar_dict, this_file) From f7e284cf88af73023e6e56116d2052b78cdd2fa7 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 11 May 2026 12:49:07 -0700 Subject: [PATCH 04/15] Update emulandice2 pipelines and fingerprintmap for Jun2024 PROTECT deliv --- modules/emulandice2/AIS/pipeline.2300.yml | 30 +++++++++++-- .../emulandice2/AIS/pipeline.global.2150.yml | 2 +- .../emulandice2/AIS/pipeline.global.2300.yml | 2 +- modules/emulandice2/AIS/pipeline.global.yml | 4 +- modules/emulandice2/GrIS/pipeline.2300.yml | 22 +++++++-- .../emulandice2/GrIS/pipeline.global.2100.yml | 2 +- .../emulandice2/GrIS/pipeline.global.2300.yml | 2 +- .../emulandice2/glaciers/pipeline.2300.yml | 6 ++- .../glaciers/pipeline.global.2300.yml | 2 +- .../emulandice2/shared/grd_fingerprintmap.yml | 45 ++++++++++++++----- 10 files changed, 90 insertions(+), 27 deletions(-) diff --git a/modules/emulandice2/AIS/pipeline.2300.yml b/modules/emulandice2/AIS/pipeline.2300.yml index 49a4cabd..42c507d8 100644 --- a/modules/emulandice2/AIS/pipeline.2300.yml +++ b/modules/emulandice2/AIS/pipeline.2300.yml @@ -27,7 +27,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_2300_240214_data.tgz" + - "emulandice2_AIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: @@ -40,7 +40,17 @@ postprocess: script_noupload: "emulandice_postprocess.py" arguments: - --ncfiles - - "%PIPELINE_ID%_ALL_globalsl.nc" + - "%PIPELINE_ID%_EAIS1_globalsl.nc" + - "%PIPELINE_ID%_EAIS2_globalsl.nc" + - "%PIPELINE_ID%_EAIS3_globalsl.nc" + - "%PIPELINE_ID%_EAIS4_globalsl.nc" + - "%PIPELINE_ID%_EAIS5_globalsl.nc" + - "%PIPELINE_ID%_EAIS6_globalsl.nc" + - "%PIPELINE_ID%_EAIS7_globalsl.nc" + - "%PIPELINE_ID%_PEN_globalsl.nc" + - "%PIPELINE_ID%_WAIS1_globalsl.nc" + - "%PIPELINE_ID%_WAIS2_globalsl.nc" + - "%PIPELINE_ID%_WAIS3_globalsl.nc" options: - "locationfile" - "chunksize" @@ -55,12 +65,24 @@ postprocess: - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' - '%MODULE_PATH%/../shared/emulandice_postprocess.py' upload_and_extract_input_data: - - 'grd_fingerprints_data.tgz' + - 'grd_fingerprints_data2.tgz' copy_shared_data: - '$SHARED/location.lst' copy_input_data: project: task1: - - "%PIPELINE_ID%_ALL_globalsl.nc" + - "%PIPELINE_ID%_EAIS1_globalsl.nc" + - "%PIPELINE_ID%_EAIS2_globalsl.nc" + - "%PIPELINE_ID%_EAIS3_globalsl.nc" + - "%PIPELINE_ID%_EAIS4_globalsl.nc" + - "%PIPELINE_ID%_EAIS5_globalsl.nc" + - "%PIPELINE_ID%_EAIS6_globalsl.nc" + - "%PIPELINE_ID%_EAIS7_globalsl.nc" + - "%PIPELINE_ID%_PEN_globalsl.nc" + - "%PIPELINE_ID%_WAIS1_globalsl.nc" + - "%PIPELINE_ID%_WAIS2_globalsl.nc" + - "%PIPELINE_ID%_WAIS3_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" local_total_files: - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/AIS/pipeline.global.2150.yml b/modules/emulandice2/AIS/pipeline.global.2150.yml index 88b52e4d..61ec5188 100644 --- a/modules/emulandice2/AIS/pipeline.global.2150.yml +++ b/modules/emulandice2/AIS/pipeline.global.2150.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_2150_240214_data.tgz" + - "emulandice2_AIS_project_2150_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/AIS/pipeline.global.2300.yml b/modules/emulandice2/AIS/pipeline.global.2300.yml index 2752d2db..6e9364f8 100644 --- a/modules/emulandice2/AIS/pipeline.global.2300.yml +++ b/modules/emulandice2/AIS/pipeline.global.2300.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_2300_240214_data.tgz" + - "emulandice2_AIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/AIS/pipeline.global.yml b/modules/emulandice2/AIS/pipeline.global.yml index af3fa527..61ec5188 100644 --- a/modules/emulandice2/AIS/pipeline.global.yml +++ b/modules/emulandice2/AIS/pipeline.global.yml @@ -10,7 +10,7 @@ project: - --region - ALL - --emu_file - - emu_file/AIS_ALL_Kori_PISM_pow_exp_10_EMULATOR.RData + - emu_file/AIS_ALL_Kori_PISM_CISM_ElmerIce_pow_exp_10_EMULATOR.RData options: - "pipeline_id" - "scenario" @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_AIS_project_data.tgz" + - "emulandice2_AIS_project_2150_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.2300.yml b/modules/emulandice2/GrIS/pipeline.2300.yml index 39d102a6..246bb416 100644 --- a/modules/emulandice2/GrIS/pipeline.2300.yml +++ b/modules/emulandice2/GrIS/pipeline.2300.yml @@ -20,6 +20,8 @@ project: - "pyear_end" - "pyear_step" - "seed" + - "cyear_start" + - "cyear_end" upload_input_data: - "%MODULE_PATH%/../shared/emulandice_environment.sh" - "%MODULE_PATH%/../shared/emulandice_steer.sh" @@ -27,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_2300_240214_data.tgz" + - "emulandice2_GrIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: @@ -40,7 +42,12 @@ postprocess: script_noupload: "emulandice_postprocess.py" arguments: - --ncfiles - - "%PIPELINE_ID%_ALL_globalsl.nc" + - "%PIPELINE_ID%_CW_globalsl.nc" + - "%PIPELINE_ID%_SW_globalsl.nc" + - "%PIPELINE_ID%_SE_globalsl.nc" + - "%PIPELINE_ID%_NE_globalsl.nc" + - "%PIPELINE_ID%_NO_globalsl.nc" + - "%PIPELINE_ID%_NW_globalsl.nc" options: - "locationfile" - "chunksize" @@ -55,12 +62,19 @@ postprocess: - '%MODULE_PATH%/../shared/grd_fingerprintmap.yml' - '%MODULE_PATH%/../shared/emulandice_postprocess.py' upload_and_extract_input_data: - - 'grd_fingerprints_data.tgz' + - 'grd_fingerprints_data2.tgz' copy_shared_data: - '$SHARED/location.lst' copy_input_data: project: task1: - - "%PIPELINE_ID%_ALL_globalsl.nc" + - "%PIPELINE_ID%_CW_globalsl.nc" + - "%PIPELINE_ID%_SW_globalsl.nc" + - "%PIPELINE_ID%_SE_globalsl.nc" + - "%PIPELINE_ID%_NE_globalsl.nc" + - "%PIPELINE_ID%_NO_globalsl.nc" + - "%PIPELINE_ID%_NW_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" local_total_files: - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/GrIS/pipeline.global.2100.yml b/modules/emulandice2/GrIS/pipeline.global.2100.yml index 1932fb9b..ef53561b 100644 --- a/modules/emulandice2/GrIS/pipeline.global.2100.yml +++ b/modules/emulandice2/GrIS/pipeline.global.2100.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_2100_240214_data.tgz" + - "emulandice2_GrIS_project_2100_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/GrIS/pipeline.global.2300.yml b/modules/emulandice2/GrIS/pipeline.global.2300.yml index 9eba19f9..8e981c11 100644 --- a/modules/emulandice2/GrIS/pipeline.global.2300.yml +++ b/modules/emulandice2/GrIS/pipeline.global.2300.yml @@ -29,7 +29,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_GrIS_project_2300_240214_data.tgz" + - "emulandice2_GrIS_project_2300_240528_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/glaciers/pipeline.2300.yml b/modules/emulandice2/glaciers/pipeline.2300.yml index d479b6b6..01e7c4a3 100644 --- a/modules/emulandice2/glaciers/pipeline.2300.yml +++ b/modules/emulandice2/glaciers/pipeline.2300.yml @@ -56,6 +56,8 @@ project: - "pyear_end" - "pyear_step" - "seed" + - "cyear_start" + - "cyear_end" upload_input_data: - "%MODULE_PATH%/../shared/emulandice_environment.sh" - "%MODULE_PATH%/../shared/emulandice_steer.sh" @@ -63,7 +65,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_glaciers_project_2300_240223_data.tgz" + - "emulandice2_glaciers_project_2300_240320_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: @@ -133,5 +135,7 @@ postprocess: - "%PIPELINE_ID%_RGI17_globalsl.nc" - "%PIPELINE_ID%_RGI18_globalsl.nc" - "%PIPELINE_ID%_RGI19_globalsl.nc" + # download_output_data: + # - "%PIPELINE_ID%_quantiles.nc" local_total_files: - "%PIPELINE_ID%_localsl.nc" diff --git a/modules/emulandice2/glaciers/pipeline.global.2300.yml b/modules/emulandice2/glaciers/pipeline.global.2300.yml index 47b0fad1..b2b657b7 100644 --- a/modules/emulandice2/glaciers/pipeline.global.2300.yml +++ b/modules/emulandice2/glaciers/pipeline.global.2300.yml @@ -65,7 +65,7 @@ project: - "%MODULE_PATH%/../shared/emulandice_project.py" upload_and_extract_input_data: - "%MODULE_PATH%/../shared/emulandice_bundled_dependencies.tgz" - - "emulandice2_glaciers_project_2300_240223_data.tgz" + - "emulandice2_glaciers_project_2300_240320_data.tgz" global_total_files: - "%PIPELINE_ID%_ALL_globalsl.nc" download_output_data: diff --git a/modules/emulandice2/shared/grd_fingerprintmap.yml b/modules/emulandice2/shared/grd_fingerprintmap.yml index a539e8a0..0b6f47da 100644 --- a/modules/emulandice2/shared/grd_fingerprintmap.yml +++ b/modules/emulandice2/shared/grd_fingerprintmap.yml @@ -6,18 +6,41 @@ # Should go to modules-data when finished AIS: - AIS.ALL: - fingerprint: FPRINT/fprint_wais.nc -# AIS.EAIS: -# fingerprint: FPRINT/fprint_eais.nc -# AIS.WAIS: -# fingerprint: FPRINT/fprint_wais.nc -# AIS.AP: -# fingerprint: FPRINT/fprint_wais.nc + AIS.WAIS1: + fingerprint: FPRINT/fprint_ais11_1.nc + AIS.WAIS2: + fingerprint: FPRINT/fprint_ais11_2.nc + AIS.WAIS3: + fingerprint: FPRINT/fprint_ais11_3.nc + AIS.EAIS1: + fingerprint: FPRINT/fprint_ais11_5.nc + AIS.EAIS2: + fingerprint: FPRINT/fprint_ais11_6.nc + AIS.EAIS3: + fingerprint: FPRINT/fprint_ais11_7.nc + AIS.PEN: + fingerprint: FPRINT/fprint_ais11_4.nc + AIS.EAIS4: + fingerprint: FPRINT/fprint_ais11_8.nc + AIS.EAIS5: + fingerprint: FPRINT/fprint_ais11_9.nc + AIS.EAIS6: + fingerprint: FPRINT/fprint_ais11_10.nc + AIS.EAIS7: + fingerprint: FPRINT/fprint_ais11_11.nc GrIS: - GIS.ALL: - fingerprint: FPRINT/fprint_gis.nc - + GIS.NW: + fingerprint: FPRINT/fprint_gis01.nc + GIS.CW: + fingerprint: FPRINT/fprint_gis02.nc + GIS.SW: + fingerprint: FPRINT/fprint_gis03.nc + GIS.SE: + fingerprint: FPRINT/fprint_gis04.nc + GIS.NE: + fingerprint: FPRINT/fprint_gis05.nc + GIS.NO: + fingerprint: FPRINT/fprint_gis06.nc glaciers: GLA.RGI01: fingerprint: FPRINT/fprint_glac11.nc From 3df3d6750ea17647fb3080ed87785d5f4d8bfeaa Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 11 May 2026 12:53:54 -0700 Subject: [PATCH 05/15] Update emulandice2 commit to May/Jun2024 --- modules/emulandice2/shared/emulandice2 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/modules/emulandice2/shared/emulandice2 b/modules/emulandice2/shared/emulandice2 index fd319046..0a646f6c 160000 --- a/modules/emulandice2/shared/emulandice2 +++ b/modules/emulandice2/shared/emulandice2 @@ -1 +1 @@ -Subproject commit fd3190469636e14bf7e6c0d1ac18a14a052fd634 +Subproject commit 0a646f6c66a86d6ffbbea9ced676236ba3405501 From b4d48ac4d2b4db734ee35ea3db8e76f2707913c1 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 11 May 2026 13:21:37 -0700 Subject: [PATCH 06/15] update protect experiment to global for now --- experiments/emulandice2.protect.ssp585/config.yml | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/experiments/emulandice2.protect.ssp585/config.yml b/experiments/emulandice2.protect.ssp585/config.yml index c1cb212f..267451e5 100644 --- a/experiments/emulandice2.protect.ssp585/config.yml +++ b/experiments/emulandice2.protect.ssp585/config.yml @@ -5,7 +5,7 @@ global-options: pyear_end: 2300 #2150 pyear_step: 10 baseyear: 2005 - pipeline_file: "pipeline.2300.yml" # "pipeline.yml" + pipeline_file: "pipeline.global.2300.yml" # "pipeline.yml" climate_step: temperature: @@ -13,6 +13,7 @@ climate_step: module: "temperature" pipeline_file: "pipeline.global.yml" generates_climate_output: true + sealevel_step: emuAIS: module_set: "emulandice2" From 3cd1cb205157d619fb91efd607aab40bf29a363f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?V=C3=ADctor=20Malag=C3=B3n=20Santos?= Date: Sat, 6 Dec 2025 01:51:02 +0100 Subject: [PATCH 07/15] EBM3 module --- .../ebm3_oceandynamics_project.py | 674 ++++++++++++++++++ .../ebm3_thermalexpansion_project.py | 164 +++++ .../sterodynamics/energy_balance_model.py | 519 ++++++++++++++ .../ebm3/sterodynamics/pipeline.global.yml | 21 + modules/ebm3/sterodynamics/pipeline.yml | 57 ++ 5 files changed, 1435 insertions(+) create mode 100644 modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py create mode 100644 modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py create mode 100644 modules/ebm3/sterodynamics/energy_balance_model.py create mode 100644 modules/ebm3/sterodynamics/pipeline.global.yml create mode 100644 modules/ebm3/sterodynamics/pipeline.yml diff --git a/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py b/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py new file mode 100644 index 00000000..e966d8c3 --- /dev/null +++ b/modules/ebm3/sterodynamics/ebm3_oceandynamics_project.py @@ -0,0 +1,674 @@ +# Sterodynamics projections for PROTECT +# Victor Malagon Santos, 17.06.2024 + +import numpy as np +import os +import pandas as pd +import sys +import argparse +import energy_balance_model as ebm3 +import csv +import xarray as xr +from sklearn import linear_model +import re +import time +import netCDF4 +import h5py +import scipy +import pickle + +""" +Created on Fri Jun 7 11:24:19 2024 + +For any model you have up tp 2300, you want a set of scaling coefficients +for ssp126 and ssp585, as those thwo a typically the ones available up to 2300. +Those scenarios are used to get a pattern of regression coefficient used for +SSP interpolation and probabilistic projections using FAIR. + +@author: vmalagonsantos +""" + +# functions +def fetch_erfs_from_rcmip(path,scenarios): + ''' + function to read scenario-dependent ERF from RCMIP csv + + Input parameters: + path = path to RCMIP csv + scenarios = desired ssp's + + Output: Adds results of extreme sea-level analysis for all stations to the + station data dictionary. + + erfs = effective radiative forcing per scenario and year + erfyears = corresponding years (1750-2500) + ''' + + erfyears=np.arange(1750,2501) + ssp_idx = {'ssp119':212,'ssp126':231,'ssp245':308,'ssp370':59,'ssp585':404} #location in table + + with open(path) as csv_file: + csv_reader = csv.reader(csv_file) + rows = list(csv_reader) + + erfs = np.empty((len(scenarios),len(erfyears))) + + for s,scen in enumerate(scenarios): + try: + erfs[s,:] = rows[ssp_idx[scen]][7:] + except: + continue + + return erfs, erfyears + +def angd(lat0, lon0, qlat, qlon): + + # Convert the input from degrees to radians + (lat0, lon0) = np.radians((lat0, lon0)) + (qlat, qlon) = np.radians((qlat, qlon)) + + # Calculate the angle between the vectors + temp = np.arctan2(np.sqrt((np.cos(qlat)*np.sin(qlon-lon0))**2 + \ + (np.cos(lat0)*np.sin(qlat) - np.sin(lat0)*np.cos(qlat) * np.cos(qlon-lon0))**2),\ + (np.sin(lat0)*np.sin(qlat) + np.cos(lat0)*np.cos(qlat)*np.cos(qlon-lon0))) + + # Convert the results from radians to degrees and return + return(np.degrees(temp)) + +def NearestPoint(qlat, qlon, lats, lons, tol = None): + + # Get the distance between the query point and all the possible points + dist = angd(lats, lons, qlat, qlon) + + # Which is the closest point + nearest_idx = np.argmin(dist) + + # Is the point within the tolerance? + if isinstance(tol, (int, float)): + if dist[nearest_idx] > tol: + return(None) + + return(nearest_idx) + + +def NearestPoints(qlats, qlons, lats, lons, tol): + + if len(qlats) != len(qlons): + raise Exception("Query lats ({}) and lons ({}) differ in length".format(len(qlats), len(qlons))) + + idx = map(lambda qlat,qlon: NearestPoint(qlat, qlon, lats, lons, tol), qlats, qlons) + + return(list(idx)) + +def ReadLocationFile(location_file): + + # Initialize variables to hold data and site information + names = [] + ids = [] + lats = [] + lons = [] + + # Compile the regex for finding commented lines + comment_regex = re.compile(r'^#') + + # Open the rate file + with open(location_file, 'r') as f: + + # Loop over the lines of the file + for line in f: + + # Skip commented lines + if re.search(comment_regex, line): + continue + + # Split the line into components + (this_name, this_id, this_lat, this_lon) = line.split("\t") + + # Store the information + names.append(this_name) + ids.append(int(this_id)) + lats.append(float(this_lat)) + lons.append(float(this_lon)) + + # Cast everything as numpy arrays + names = np.array(names) + ids = np.array(ids) + lats = np.array(lats) + lons = np.array(lons) + + # Return variables + return(names, ids, lats, lons) + + +# Smooth ZOS and ZOSTOGA over 19 year smoothing window +def nanSmooth(x, w=19): + idx = np.flatnonzero(~np.isnan(x)) + temp = x + if len(idx) > 0: + temp[idx] = Smooth(x[idx], w) + return(temp) + +def Smooth(x, w=19): + out0 = np.convolve(x, np.ones(w,dtype='double'), 'valid')/w + r = np.arange(1,w-1,2, dtype="double") + start = np.cumsum(x[:w-1])[::2]/r + stop = (np.cumsum(x[:-w:-1])[::2]/r)[::-1] + y = np.concatenate((start, out0, stop)) + return(y) + +def emb3_thermalexpansion_postprocess(scenario, pipeline_id, nsamps, seed, pyear_start, pyear_end, pyear_step, locationfile, baseyear, climate_data_file, rfmip, params, zosdir): + if pyear_end < 2151: + ryear_end = 2100 + elif pyear_end > 2300: + ryear_end = 2300 + else: + ryear_end = pyear_end + + targyears = np.arange(pyear_start, ryear_end) # for regression + projyears = np.arange(pyear_start, pyear_end+1, pyear_step) + + (_, site_ids, site_lats, site_lons) = ReadLocationFile(locationfile) + + # get temperature from FaIR simulations: INPUT FROM CLIMATE STEP. We need both gsat and oceantemp + cf = xr.open_dataset(climate_data_file, group=scenario, engine='netcdf4', ) + gsat = cf['surface_temperature'].sel(years=projyears) - cf['surface_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + otemp = cf['deep_ocean_temperature'].sel(years=projyears) - cf['deep_ocean_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + # INPUT temperature file from ebm3 global + gte_file = f'{pipeline_id}_globalsl.nc' + + # INPUT get models and parameters + paramdir = params + + #forcing for ebm + scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections + rfmipfile = rfmip #path # RFMIP FILE IS INPUT, NEW DATA + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries + + ebm_param = pd.read_csv(paramdir) + + # find available zos simulations - ignoring hidden files + model_has_zos = [f for f in os.listdir(zosdir) if not f.startswith('.')] + + ## FITTING + slopes = [] + intercept = [] + + # if projections are below 2150, look for 2100 zos simulations available for requested SSP, run ebm adn estimate slopes + if pyear_end < 2151: # # keep those up to 2300 or 2500 + for m, model in enumerate(model_has_zos): + ssp1_file = [] + hist_file = [] + # model = 'ACCESS-ESM1-5' + + # find matching ebm parameters to model + # df['company_name'].eq('ABC').any() + model_param = ebm_param[ebm_param["model"].eq(model)] + modeldir = zosdir + model + + zos_runs2_2100 = [f for f in os.listdir(modeldir) if scenario in f] # regex would be better here + if not zos_runs2_2100: + print(' ') + print(scenario + ' scenario not found in model ' + model + '. Moving on.') + continue + print(' ') + print('Estimating regression parameters for model: ' + model) + + + # get variants for 2100 simulation + + for v, filename in enumerate(zos_runs2_2100): + variant = [f for f in re.split('_', filename) if f.startswith('r')][0] + # get EBM parameters for both model and variant + model_var_param = model_param[model_param["run"].str.contains(variant)] + if not model_var_param.empty: #if not empty, check if that variant is availble for both hist and ssps + # get zos, historical + ssp + hist_file = [f for f in os.listdir(modeldir) if 'hist' in f and variant in f] + ssp_file = [f for f in os.listdir(modeldir) if scenario in f and variant in f] + if (hist_file and ssp_file): # if there is zos sims for hist and ssps and ebm parameters, continue with emulation + break + if model_var_param.empty: + print('No EBM parameters found for model ' + model) + continue + + + # get model parameters + c1 = np.array(model_var_param.get('C1')) + c2 = np.array(model_var_param.get('C2')) + c3 = np.array(model_var_param.get('C3')) + k1 = np.array(model_var_param.get('kappa1')) + k2 = np.array(model_var_param.get('kappa2')) + k3 = np.array(model_var_param.get('kappa3')) + e = np.array(model_var_param.get('epsilon')) + f4 = np.array(model_var_param.get('F_4xCO2')) + se = np.array(model_var_param.get('sigma_eta')) + sx = np.array(model_var_param.get('sigma_xi')) + gamma = np.array(model_var_param.get('gamma')) + nit = np.array(model_var_param.get('nit')) + + + # get forcing + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, [scenario]) #get ERF timeseries + + + #Model setup 3lm + ebm_obj = ebm3.EnergyBalanceModel( + ocean_heat_capacity = [c1[0], c2[0], c3[0]], + ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value + deep_ocean_efficacy = e[0], + forcing_4co2 = f4[0], + stochastic_run=False, + sigma_eta = se[0], + sigma_xi = sx[0], + gamma_autocorrelation = gamma[0], + seed = None, + timestep = 1, + n_timesteps = 2500-1750, + ) + + # Run model for different scenarios + ebm_obj.add_forcing( erfs[0,:], erfyears) + ebm_obj.run() + temp3 = ebm_obj.temperature + + attrs={'description': 'surface and deep temperature obtained by EBM-3LM', + 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} + + temp3_xr = xr.Dataset({'temperature':(['years','layers'], temp3)}, + coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears}, attrs=attrs) + + # reference to base year + temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) + # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + try: + hist_xr = xr.open_dataset(modeldir + '/' + hist_file[0],decode_times=False) + ssp_xr = xr.open_dataset(modeldir + '/' + ssp_file[0],decode_times=False) + except: + print('Model misses either historical or ssp simulation. Moving on') + continue + + if ssp_xr['time'][0].values != 60280.5: + ssp_xr.coords['time'] = ((ssp_xr['time'] + 60280.5).astype('float')) + + # cocatenate + zos = xr.concat([hist_xr, ssp_xr], dim='time') + + #% create new coordinates for year + zos.coords['years'] = ((zos['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + + # Calculate annual means + zos_am = zos.groupby('years').mean(dim='time') + zos_am = zos_am.where(zos_am['zos'] < 99999) # replace land values with nans + + # reference zos + # baseyear_idx = np.flatnonzero(datayears == baseyear) + # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) + try: + zos_am = zos_am.sel(years=targyears) - zos_am.sel(years=baseyear) + except: + print('Model has inconsistent time. Moving on') + continue + + #% fitting stage + nlon = len(zos_am['lon']) + nlat = len(zos_am['lat']) + + # get temperature for desired scenario + Ts = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='surface')) + Ti = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='inter')) + Td = np.array(temp3_xr['temperature'].sel(years=targyears).sel(layers='deep')) + + + zos = np.array(zos_am['zos'].sel(years=targyears)) + + slope = np.zeros((3, nlat, nlon)) + slope[:] = np.nan + smoothwin = 19 + + reg = linear_model.LinearRegression() + # for i in tqdm(range(nlat)): + # sleep(3) + for i in range(nlat): + for j in range(nlon): + + y = np.array(Smooth(zos[:,i,j],w=smoothwin)) # reducing varability + x1 = Ts.flatten() + x2 = Ti.flatten() + x3 = Td.flatten() + x = np.transpose((x1,x2,x3)) + + try: + + reg.fit(x, y) + slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface + slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter + slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep + # intercept[m,i,j] = reg.intercept_ + + except ValueError: + slope[0,i,j] = np.nan + slope[1,i,j] = np.nan + slope[2,i,j] = np.nan + # intercept[m,i,j] = np.nan + + slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? + + + # if projections are above 2150, look for 2300 zos simulations available for SSP126 and SSP585, run ebm, estimate slopes, and emulate requested scenario + if pyear_end >= 2151: + for m, model in enumerate(model_has_zos): + + # find matching ebm parameters to model + model_param = ebm_param[ebm_param["model"].str.contains(model)] + modeldir = zosdir + model + + # keep those up to 2300 or 2500 + zos_runs2_2300 = [f for f in os.listdir(modeldir) if '2300' in f or '2500' in f] # regex would be better here + if not zos_runs2_2300: + print(' ') + print('Model ' + model + ' only runs to 2100. Skipping and looking for longer simulations' ) + continue + print(' ') + print('Estimating parameters for model: ' + model) + + for v, filename in enumerate(zos_runs2_2300): + ssp1_file = [] + ssp5_file = [] + hist_file = [] + + variant = [f for f in re.split('_', filename) if f.startswith('r')][0] + + # get EBM parameters for both model and variant + model_var_param = model_param[model_param["run"].str.contains(variant)] + if not model_var_param.empty: #if not empty, check if that variant is availble for both hist and ssps + # get zos, historical + ssp + hist_file = [f for f in os.listdir(modeldir) if 'hist' in f and variant in f] + ssp1_file = [f for f in os.listdir(modeldir) if 'ssp126' in f and variant in f] + ssp5_file = [f for f in os.listdir(modeldir) if 'ssp585' in f and variant in f] + if (hist_file and ssp1_file and ssp5_file): # if there is zos sims for hist and ssps and ebm parameters, continue with emulation + break + + + # define EBM parameters + c1 = np.array(model_var_param.get('C1')) + c2 = np.array(model_var_param.get('C2')) + c3 = np.array(model_var_param.get('C3')) + k1 = np.array(model_var_param.get('kappa1')) + k2 = np.array(model_var_param.get('kappa2')) + k3 = np.array(model_var_param.get('kappa3')) + e = np.array(model_var_param.get('epsilon')) + f4 = np.array(model_var_param.get('F_4xCO2')) + se = np.array(model_var_param.get('sigma_eta')) + sx = np.array(model_var_param.get('sigma_xi')) + gamma = np.array(model_var_param.get('gamma')) + nit = np.array(model_var_param.get('nit')) + + scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections + erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries + + + temp3 = np.empty((len(erfyears), 3, len(scenarios))) + + for s, scen in enumerate(scenarios): + + #Model setup 3lm + ebm_obj = ebm3.EnergyBalanceModel( + ocean_heat_capacity = [c1[0], c2[0], c3[0]], + ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value + deep_ocean_efficacy = e[0], + forcing_4co2 = f4[0], + stochastic_run=False, + sigma_eta = se[0], + sigma_xi = sx[0], + gamma_autocorrelation = gamma[0], + seed = None, + timestep = 1, + n_timesteps = 2500-1750, + ) + + # Run model for different scenarios + ebm_obj.add_forcing( erfs[s,:], erfyears) + ebm_obj.run() + temp3[:,:,s] = ebm_obj.temperature + + attrs={'description': 'surface and deep temperature obtained by EBM-3LM', + 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} + + temp3_xr = xr.Dataset({'temperature':(['years','layers', 'scenarios'], temp3)}, + coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears, 'scenarios': scenarios}, attrs=attrs) + + # reference to base year + temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) + # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') + + try: + hist_xr = xr.open_dataset(modeldir + '/' + hist_file[0],decode_times=False) + ssp1_xr = xr.open_dataset(modeldir + '/' + ssp1_file[0],decode_times=False) + ssp5_xr = xr.open_dataset(modeldir + '/' + ssp5_file[0],decode_times=False) + except: + print('Model misses one of the scenarios used for emulation. Moving on.') + continue + + # some models' time is references to 2015 instead of 1850 + if ssp1_xr['time'][0].values != 60280.5: + ssp1_xr.coords['time'] = ((ssp1_xr['time'] + 60280.5).astype('float')) + ssp5_xr.coords['time'] = ((ssp5_xr['time'] + 60280.5).astype('float')) + + # cocatenate + zos1 = xr.concat([hist_xr, ssp1_xr], dim='time') + zos5 = xr.concat([hist_xr, ssp5_xr], dim='time') + + #% create new coordinates for year + zos1.coords['years'] = ((zos1['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + zos5.coords['years'] = ((zos5['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR + + # Calculate annual means + zos1_am = zos1.groupby('years').mean(dim='time') + zos1_am = zos1_am.where(zos1_am['zos'] < 99999) # replace land values with nans + zos5_am = zos5.groupby('years').mean(dim='time') + zos5_am = zos5_am.where(zos5_am['zos'] < 99999) # replace land values with nans + + # reference zos + # baseyear_idx = np.flatnonzero(datayears == baseyear) + # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) + zos1_am = zos1_am.sel(years=targyears) - zos1_am.sel(years=baseyear) + zos5_am = zos5_am.sel(years=targyears) - zos5_am.sel(years=baseyear) + + #% fitting stage + nlon = len(zos1_am['lon']) + nlat = len(zos1_am['lat']) + + + Ts = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='surface').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='surface').drop_vars('scenarios')),dim='years') + + Ti = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='inter').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='inter').drop_vars('scenarios')),dim='years') + + Td = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='deep').drop_vars('scenarios'), + temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='deep').drop_vars('scenarios')),dim='years') + + + zos1 = np.array(zos1_am['zos'].sel(years=targyears)) + zos5 = np.array(zos5_am['zos'].sel(years=targyears)) + + slope = np.zeros((3, nlat, nlon)) + slope[:] = np.nan + smoothwin = 19 + + reg = linear_model.LinearRegression() + # for i in tqdm(range(nlat)): + # sleep(3) + for i in range(nlat): + for j in range(nlon): + + y = np.concatenate((Smooth(zos1[:,i,j],w=smoothwin), Smooth(zos5[:,i,j],w=smoothwin))) # reducing varability + x1 = np.array(Ts['temperature'][:]).flatten() + x2 = np.array(Ti['temperature'][:]).flatten() + x3 = np.array(Td['temperature'][:]).flatten() + x = np.transpose(np.stack((x1,x2,x3))) + + try: + + reg.fit(x, y) + slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface + slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter + slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep + # intercept[m,i,j] = reg.intercept_ + + except ValueError: + slope[0,i,j] = np.nan + slope[1,i,j] = np.nan + slope[2,i,j] = np.nan + # intercept[m,i,j] = np.nan + + slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? + + lons, lats = np.meshgrid(hist_xr['lon'], hist_xr['lat']) + lats = lats.flatten() + lons = lons.flatten() + + site_ids_map = np.array(NearestPoints(site_lats, site_lons, lats, lons, tol=None)) + + # # remove outliers + # slopes_sum = np.empty((len(slopes), len(lats))) + # # slopes_sum = np.array([np.sum(s,0).flatten() for s in slopes]) + # slopes_sum = np.array([s[0,:,:].flatten() for s in slopes]) + # slopes_sum_std = np.std(slopes_sum,0) + # [ev, ef] = np.unique(np.where(slopes_sum > slopes_sum_std * 3 )[0], return_counts = True) + # mi = [i for i,v in enumerate(ef) if v > 1000] + # slopes = [sl for s, sl in enumerate(slopes) if not s in mi ] + + #nsamps + samples = np.array(gsat['samples']) + + # RESAMPLE SLOPE PARAMETERS + nsims=len(slopes) + + rng = np.random.default_rng(seed) + if nsamps > nsims: + run_idx = np.arange(nsims) + sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) + else: + run_idx = rng.choice(nsims, nsamps, nsamps>nsims) + sample_idx = np.arange(nsamps) + + + ## resample slopes + slopes_resampled = [slopes[i] for i in sample_idx] + + dsl = np.empty((nsamps, len(projyears), len(site_ids))) + dsl[:] = np.nan + + for sample in samples: + slope = slopes_resampled[sample] + print(sample) + + # prepare slopes + slope_s = slope[0,:,:].flatten() + slope_s = np.array([slope_s[x] for x in site_ids_map]) + + slope_i = slope[1,:,:].flatten() + slope_i = np.array([slope_i[x] for x in site_ids_map]) + + slope_d = slope[2,:,:].flatten() + slope_d = np.array([slope_d[x] for x in site_ids_map]) + + # obtain fair temepratures + Tfs = np.array(gsat.sel(samples=sample)) # surface temperature from fair + Tfi = np.array(otemp.sel(layers=1).sel(samples=sample)) # intermediate temperature from fair + Tfd = np.array(otemp.sel(layers=2).sel(samples=sample)) # deep temperature from fair + + + #% project + dsl[sample,:,:] = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) + + # dsl = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) + + ncvar_attributes = {"description": "Dynamic Sea Level"} + + nc_missing_value = np.nan + # Generate the output xarray + dsl_xr = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), dsl, {"units":"mm", "missing_value":nc_missing_value}), + "lat": (("locations"), site_lats), + "lon": (("locations"), site_lons)}, + coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) + + + # make sure projections have 0 mean. This weighted average only works for regular grids. + # Must be edited to accomodate other grids if zos is not in 1x1. + + weights = np.cos(np.deg2rad(dsl_xr.lat)) + weights.name = "weights" + + dsl_xr_weighted = dsl_xr.weighted(weights) + dsl_xr_weighted + weighted_mean = dsl_xr_weighted.mean(("locations")) + + + # Appply 0mean correction to dsl, and add GTE to get sterodynamics + + gte = xr.open_dataset(gte_file) + + + gte = gte['sea_level_change'].values[:,:,0] + sdsl = dsl*1000 - np.repeat(np.array(weighted_mean['sea_level_change'])[:, :, np.newaxis], len(site_ids), axis=2) + np.repeat(gte[0:nsamps, :, np.newaxis], len(site_ids), axis=2) + + ncvar_attributes = {"description": "Local SLR contributions from thermal expansion and dynamic sea-level using EBM3", + "history": "Created " + time.ctime(time.time()), + "source": "SLR Framework: PROTECT 2300", + "scenario": scenario, # change to scenario later + "baseyear": baseyear} + + nc_missing_value = np.nan + # Generate the output xarray + local_out = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), sdsl, {"units":"mm", "missing_value":nc_missing_value}), + "lat": (("locations"), site_lats), + "lon": (("locations"), site_lons)}, + coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) + # Write these samples to a temporary netcdf file + local_out.to_netcdf("{0}_localsl.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) + local_outq = local_out.quantile([0.01,0.05,0.17,0.50,0.83,0.95,0.99], dim='samples') + local_outq.to_netcdf("{0}_quantiles.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) + +if __name__ == '__main__': + + # Initialize the command-line argument parser + parser = argparse.ArgumentParser(description="Compute dynamic sea level for EBM3 workflow",\ + epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") + + # Define the command line arguments to be expected + parser.add_argument('--scenario', help="SSP scenario (i.e ssp585)", default='ssp585') + parser.add_argument('--nsamps', help="Number of samples to generate [default=20000]", default=20000, type=int) + parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) + parser.add_argument('--pyear_start', help="Year for which projections start [default=2000]", default=2000, type=int) + parser.add_argument('--pyear_end', help="Year for which projections end [default=2300]", default=2300, type=int) + parser.add_argument('--pyear_step', help="Step size in years between pyear_start and pyear_end at which projections are produced [default=5]", default=5, type=int) + parser.add_argument('--locationfile', help="File that contains name, id, lat, and lon of points for localization", default="location.lst") + parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) + parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") + parser.add_argument('--climate_data_file',type=str) + parser.add_argument('--rfmip', help='rfmip file',default='rfmip-radiative-forcing-annual-means-v4-0-0.csv') + parser.add_argument('--params', help='CMIP6 Params cvs', default='4xCO2_cummins_ebm3_cmip6.csv') + parser.add_argument('--zosdir',help='Path to CMIP6 ZOS directory', default='cmip6/zos/') + + # Parse the arguments + args = parser.parse_args() + + emb3_thermalexpansion_postprocess(args.scenario, + args.pipeline_id, + args.nsamps, + args.seed, + args.pyear_start, + args.pyear_end, + args.pyear_step, + args.locationfile, + args.baseyear, + args.climate_data_file, + args.rfmip, + args.params, + args.zosdir) + + # Done + sys.exit() \ No newline at end of file diff --git a/modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py b/modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py new file mode 100644 index 00000000..564ce43a --- /dev/null +++ b/modules/ebm3/sterodynamics/ebm3_thermalexpansion_project.py @@ -0,0 +1,164 @@ +# Thermal expansion projections for PROTECT +# Victor Malagon Santos, 17.06.2024 + +import os +import numpy as np +import netCDF4 as nc +import pandas as pd +import time +from netCDF4 import Dataset +import pickle +import sys +import re +import argparse + +class ProjectionError(Exception): + pass + + +def ebm3_thermalexpansion_project(scenario, climate_data_file, coef_file, params_file, pyear_start, pyear_end, pyear_step, nsamps, pipeline_id, seed, baseyear): + # constants + a = 6.37*1e6 + earth_area = 4*np.pi*a**2 + + path = os.path.dirname(__file__) + + # heat capacity of each layer from fair2. Needed to compute OHC + fparam = pd.read_csv(params_file) + c1 = fparam['clim_c1'] + c2 = fparam['clim_c2'] + c3 = fparam['clim_c3'] + + nsims=1001 + + rng = np.random.default_rng(seed) + if nsamps > nsims: + run_idx = np.arange(nsims) + sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) + else: + run_idx = rng.choice(nsims, nsamps, nsamps>nsims) + sample_idx = np.arange(nsamps) + + c1 = c1[sample_idx] + c2 = c2[sample_idx] + c3 = c3[sample_idx] + + # create target years array + targyears = np.arange(pyear_start, pyear_end+1, pyear_step) + + #temperature output from fair2. Needed t compute OHC + ds = Dataset(climate_data_file) + gsat = ds[scenario]['surface_temperature'][:] + deepoceant = ds[scenario]['deep_ocean_temperature'][:] + years = ds[scenario]['years'][:] + + # Expansion coefficients. Needed to to convert OCH to global thermal expansion (GTE) + ds = Dataset(coef_file) + include_models = ds['model'][:] + eeh3 = ds['expcoefs'][:] + + # estimate OHC + up = gsat*np.array(c1) + mid = deepoceant[:,:,0]*np.array(c2) + low = deepoceant[:,:,1]*np.array(c3) + n = np.zeros((len(years), nsamps)) + n[1:,:] = up[1:,:]-up[:-1,:] + mid[1:,:]-mid[:-1,:] + low[1:,:]-low[:-1,:] + ohc_samps = np.cumsum(n*earth_area, 0) * 365*24*3.6*1e3*1e-24 + + # Generate samples assuming normal distribution + rng = np.random.default_rng(seed) + expcoef_samps = rng.normal(loc=np.mean(eeh3), scale=np.std(eeh3), size=(nsamps,1)) + + # Produce the projection samples + gte_samps = ohc_samps * expcoef_samps.flatten() + + # Center these samples on the baseyear + baseyear_idx = np.flatnonzero(years == baseyear) + gte_samps = gte_samps - gte_samps[baseyear_idx,:] + + # Subset the samples for the projection years + targyear_idx = np.isin(years, targyears) + gte_samps = gte_samps[targyear_idx,:] + + # Invert the dimensions of the variable and convert from m to mm + gte_samps *= 1000. + gte_samps = gte_samps.T + + # Save the projections to a pickle + output = {"thermsamps": gte_samps, "targyears": targyears, "baseyear": baseyear, \ + "include_models": include_models, "scenario": scenario} + outfile = open(os.path.join(os.path.dirname(__file__), "{}_projections.pkl".format(pipeline_id)), 'wb') + pickle.dump(output, outfile) + outfile.close() + + # Write the total global projections to a netcdf file + nc_filename = os.path.join(os.path.dirname(__file__), "{0}_globalsl.nc".format(pipeline_id)) + rootgrp = Dataset(nc_filename, "w", format="NETCDF4") + + # Define Dimensions + nyr = len(targyears) + year_dim = rootgrp.createDimension("years", nyr) + samp_dim = rootgrp.createDimension("samples", nsamps) + loc_dim = rootgrp.createDimension("locations", 1) + + # Populate dimension variables + year_var = rootgrp.createVariable("years", "i4", ("years",)) + samp_var = rootgrp.createVariable("samples", "i8", ("samples",)) + loc_var = rootgrp.createVariable("locations", "i8", ("locations",)) + lat_var = rootgrp.createVariable("lat", "f4", ("locations",)) + lon_var = rootgrp.createVariable("lon", "f4", ("locations",)) + + # Create a data variable + samps = rootgrp.createVariable("sea_level_change", "f4", ("samples", "years", "locations"), zlib=True, complevel=4) + + # Assign attributes + rootgrp.description = "Global SLR contribution from Thermal Expansion according to Three-Layer Model workflow" + rootgrp.history = "Created " + time.ctime(time.time()) + rootgrp.source = "FACTS: {0}".format(pipeline_id) + rootgrp.scenario = scenario + rootgrp.baseyear = baseyear + rootgrp.comment = "Included Models: " + ",".join([str(x) for x in include_models]) + samps.units = "mm" + + # Put the data into the netcdf variables + year_var[:] = targyears + samp_var[:] = np.arange(nsamps) + samps[:,:,:] = gte_samps[:,:,np.newaxis] + lat_var[:] = np.inf + lon_var[:] = np.inf + loc_var[:] = -1 + + # Close the netcdf + rootgrp.close() + + return(0) + + + +if __name__ == '__main__': + + # Initialize the command-line argument parser + parser = argparse.ArgumentParser(description="Run the global mean thermosteric sea-level rise projections",\ + epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") + + # Define the command line arguments to be expected + parser.add_argument('--scenario', help="SSP scenario (i.e ssp585) or temperature target (i.e. tlim2.0win0.25)", default='ssp585') + parser.add_argument('--climate_data_file', help="NetCDF4/HDF5 file containing surface temperature data", type=str) + parser.add_argument('--params_file', help='Full path to calibrated constraints params file', default='calibrated_constrained_parameters.csv') + parser.add_argument('--coef_file', help='Full path to expansion coefficient file', default='scmpy3LM_RCMIP_CMIP6calpm_n18_expcoefs.nc') + parser.add_argument('--nsamps', help="Number of samples to generate [default=1000]", default=1000, type=int) + parser.add_argument('--pyear_start', help="Projection year start [default=2020]", default=2020, type=int) + parser.add_argument('--pyear_end', help="Projection year end [default=2100]", default=2150, type=int) + parser.add_argument('--pyear_step', help="Projection year step [default=10]", default=10, type=int) + parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) + parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) + parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") + + + # Parse the arguments + args = parser.parse_args() + + # Run the projection process on the files specified from the command line argument + ebm3_thermalexpansion_project(args.scenario, args.climate_data_file, args.coef_file, args.params_file, args.pyear_start, args.pyear_end, args.pyear_step, args.nsamps, args.pipeline_id, args.seed, args.baseyear) + + exit() diff --git a/modules/ebm3/sterodynamics/energy_balance_model.py b/modules/ebm3/sterodynamics/energy_balance_model.py new file mode 100644 index 00000000..b110180b --- /dev/null +++ b/modules/ebm3/sterodynamics/energy_balance_model.py @@ -0,0 +1,519 @@ +"""n-layer energy balance representation of Earth's climate.""" + +import numpy as np +import scipy.linalg +import scipy.sparse.linalg +import scipy.stats +import xarray as xr + +DOUBLING_TIME_1PCT = np.log(2) / np.log(1.01) + +#: Radius of Earth, m +earth_radius = 6371000 + +#: Mass of Earth's atmosphere, kg. +mass_atmosphere = 5.1352e18 + +#: Length of a tropical year, s. +seconds_per_year = 60 * 60 * 24 * 365.24219 + + +class EnergyBalanceModel: + """Energy balance model that converts forcing to temperature. + + The energy balance model is converted to an impulse-response formulation + (hence the IR part of FaIR) to allow efficient evaluation. The benefits of + this are increased as once derived, the "layers" of the energy balance + model do not communicate with each other. The model description can be + found in [Leach2021]_, [Cummins2020]_, [Tsutsui2017]_ and [Geoffroy2013]_. + + Parameters + ---------- + ocean_heat_capacity : ``np.ndarray`` + Ocean heat capacity of each layer (top first), W m-2 yr K-1 + ocean_heat_transfer : ``np.ndarray`` + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : float + efficacy of deepest ocean layer. See e.g. [Geoffroy2013]_. + forcing_4co2 : float + effective radiative forcing from a quadrupling of atmospheric + CO2 concentrations above pre-industrial. + stochastic_run : bool + Activate the stochastic variability component from [Cummins2020]_. + sigma_eta : float + Standard deviation of stochastic forcing component from [Cummins2020]_. + sigma_xi : float + Standard deviation of stochastic disturbance applied to surface + layer. See [Cummins2020]_. + gamma_autocorrelation : float + Stochastic forcing continuous-time autocorrelation parameter. + See [Cummins2020]_. + seed : int or None + Random seed to use for stochastic variability. + timestep : float + Time interval of the model (yr) + + Raises + ------ + ValueError + if the shapes of ``ocean_heat_capacity`` and ``ocean_heat_transfer`` differ. + ValueError + if there are not at least two layers in the energy balance model. + """ + + def __init__( + self, + ocean_heat_capacity, + ocean_heat_transfer, + deep_ocean_efficacy=1, + forcing_4co2=8, + stochastic_run=False, + sigma_eta=0.5, + sigma_xi=0.5, + gamma_autocorrelation=2, + seed=None, + timestep=1, + n_timesteps=1, + ): + """Initialise the EnergyBalanceModel.""" + # adjust ocean heat capacity to be a rate: units W m-2 K-1 + self.ocean_heat_transfer = np.asarray(ocean_heat_transfer) + self.deep_ocean_efficacy = deep_ocean_efficacy + self.forcing_4co2 = forcing_4co2 + self.stochastic_run = stochastic_run + self.sigma_eta = sigma_eta + self.sigma_xi = sigma_xi + self.gamma_autocorrelation = gamma_autocorrelation + self.ocean_heat_capacity = np.asarray(ocean_heat_capacity) / timestep + self.n_temperature_boxes = len(self.ocean_heat_capacity) + if len(self.ocean_heat_transfer) != self.n_temperature_boxes: + raise ValueError( + "ocean_heat_capacity and ocean_heat_transfer must be arrays of the " + "same shape." + ) + if self.n_temperature_boxes < 2: + raise ValueError( + "At least two ocean layers must be specified in the energy balance " + "model." + ) + self.temperature = np.zeros((1, self.n_temperature_boxes + 1)) + self.n_timesteps = n_timesteps + self.n_matrix = self.n_temperature_boxes + 1 + self.seed = seed + self.timestep = timestep + + def _eb_matrix(self): + """Define the matrix of differential equations. + + Returns + ------- + eb_matrix_eigenvalues : `np.ndarray` + 1D array of eigenvalues of the energy balance matrix. + eb_matrix_eigenvectors : `np.ndarray` + 2D array of eigenvectors (an array of 1D eigenvectors) of the + energy balance matrix. + """ + # two box model + # [x x] + # [x x] + + # three box model + # [x x 0] + # [x x ex] + # [0 x x] + + # four box model + # [x x 0 0] + # [x x x 0] + # [0 x x ex] + # [0 0 x x] + + # put the efficacy of deep ocean in the right place + # making a vector avoids if statements + n_box = self.n_temperature_boxes + eb_matrix = np.zeros((n_box, n_box)) + epsilon_array = np.ones(n_box) + epsilon_array[n_box - 2] = self.deep_ocean_efficacy + + # First row + eb_matrix[0, :2] = [ + -( + self.ocean_heat_transfer[0] + + epsilon_array[0] * self.ocean_heat_transfer[1] + ) + / self.ocean_heat_capacity[0], + epsilon_array[0] + * self.ocean_heat_transfer[1] + / self.ocean_heat_capacity[0], + ] + + # Last row + eb_matrix[-1, -2:] = [ + self.ocean_heat_transfer[-1] / self.ocean_heat_capacity[-1], + -self.ocean_heat_transfer[-1] / self.ocean_heat_capacity[-1], + ] + + # Intermediate rows where n>2 + for row in range(1, n_box - 1): + eb_matrix[row, row - 1 : row + 2] = [ + self.ocean_heat_transfer[row] / self.ocean_heat_capacity[row], + -( + self.ocean_heat_transfer[row] + + epsilon_array[row] * self.ocean_heat_transfer[row + 1] + ) + / self.ocean_heat_capacity[row], + epsilon_array[row] + * self.ocean_heat_transfer[row + 1] + / self.ocean_heat_capacity[row], + ] + + # Prepend eb_matrix with stochastic terms if this is a stochastic run: + # Cummins et al. (2020) eqs. 13 and 14 + eb_matrix = np.insert(eb_matrix, 0, np.zeros(n_box), axis=0) + prepend_col = np.zeros(n_box + 1) + prepend_col[0] = -self.gamma_autocorrelation + prepend_col[1] = 1 / self.ocean_heat_capacity[0] + eb_matrix = np.insert(eb_matrix, 0, prepend_col, axis=1) + return eb_matrix + + @property + def eb_matrix_d(self): + """Return the discretised matrix exponential.""" + _eb_matrix_d = scipy.linalg.expm(self._eb_matrix()) + return _eb_matrix_d + + def _forcing_vector(self): + forcing_vector = np.zeros(self.n_temperature_boxes + 1) + forcing_vector[0] = self.gamma_autocorrelation + return forcing_vector + + @property + def forcing_vector_d(self): + """Return the discretised forcing vector.""" + return scipy.linalg.solve( + self._eb_matrix(), + (self.eb_matrix_d - np.identity(self.n_temperature_boxes + 1)) + @ self._forcing_vector(), + ) + + @property + def stochastic_d(self): + """Return the stochastic matrix.""" + # define stochastic matrix + _stochastic_d = np.zeros((self.n_timesteps, self.n_temperature_boxes + 1)) + + # stochastic stuff + if self.stochastic_run: + eb_matrix = self._eb_matrix() + q_mat = np.zeros((self.n_matrix, self.n_matrix)) + q_mat[0, 0] = self.sigma_eta**2 + q_mat[1, 1] = (self.sigma_xi / self.ocean_heat_capacity[0]) ** 2 + # use Van Loan (1978) to compute the matrix exponential + h_mat = np.zeros((self.n_matrix * 2, self.n_matrix * 2)) + h_mat[: self.n_matrix, : self.n_matrix] = -eb_matrix + h_mat[: self.n_matrix, self.n_matrix :] = q_mat + h_mat[self.n_matrix :, self.n_matrix :] = eb_matrix.T + g_mat = scipy.sparse.linalg.expm(h_mat) + q_mat_d = ( + g_mat[self.n_matrix :, self.n_matrix :].T + @ g_mat[: self.n_matrix, self.n_matrix :] + ) + q_mat_d = q_mat_d.astype(np.float64) + _stochastic_d = scipy.stats.multivariate_normal.rvs( + size=self.n_timesteps, + mean=np.zeros(self.n_matrix), + cov=q_mat_d, + random_state=self.seed, + ) + + return _stochastic_d + + def impulse_response(self): + """Convert the energy balance to impulse response representation.""" + eb_matrix = self._eb_matrix() + + # calculate the eigenvectors and eigenvalues on the energy balance + # (determininstic) part of the matrix, these are the timescales of responses + eb_matrix_eigenvalues, eb_matrix_eigenvectors = scipy.linalg.eig( + eb_matrix[1:, 1:] + ) + self.timescales = -self.timestep / (np.real(eb_matrix_eigenvalues)) + self.response_coefficients = ( + self.timescales + * ( + eb_matrix_eigenvectors[0, :] + * scipy.linalg.inv(eb_matrix_eigenvectors)[:, 0] + ) + / (self.ocean_heat_capacity[0] * self.timestep) + ) + + def emergent_parameters(self, forcing_2co2_4co2_ratio=0.5): + """Calculate emergent parameters from the energy balance parameters. + + Parameters + ---------- + forcing_2co2_4co2_ratio : float + ratio of (effective) radiative forcing converting a quadrupling of + CO2 to a doubling of CO2. + """ + # requires impulse response step + if not hasattr(self, "timescales"): + self.impulse_response() + self.ecs = ( + self.forcing_4co2 + * forcing_2co2_4co2_ratio + * np.sum(self.response_coefficients) + ) + self.tcr = ( + self.forcing_4co2 + * forcing_2co2_4co2_ratio + * np.sum( + self.response_coefficients + * ( + 1 + - self.timescales + / DOUBLING_TIME_1PCT + * (1 - np.exp(-DOUBLING_TIME_1PCT / self.timescales)) + ) + ) + ) + + def add_forcing(self, forcing, timestep): + """Add a forcing time series to EnergyBalanceModel. + + Parameters + ---------- + forcing : np.ndarray + time series of [effective] radiative forcing + timestep : float + Model timestep, in years + """ + self.forcing = forcing + self.timestep = timestep + self.n_timesteps = len(forcing) + + def run(self): + """Run the EnergyBalanceModel.""" + # internal variables + n_box = self.n_matrix - 1 + forcing_vector = self._forcing_vector() + + # Calculate the matrix exponential + eb_matrix = self._eb_matrix() + eb_matrix_d = scipy.linalg.expm(eb_matrix) + + # Solve for temperature + forcing_vector_d = scipy.linalg.solve( + eb_matrix, (eb_matrix_d - np.identity(self.n_matrix)) @ forcing_vector + ) + + solution = np.zeros((self.n_timesteps, self.n_matrix)) + solution[0, :] = self.temperature[0, :] + for i in range(1, self.n_timesteps): + solution[i, :] = ( + eb_matrix_d @ solution[i - 1, :] + + forcing_vector_d * self.forcing[i - 1] + + self.stochastic_d[i - 1, :] + ) + + self.temperature = solution[:, 1:] + self.stochastic_forcing = solution[:, 0] + self.toa_imbalance = ( + self.forcing + - self.ocean_heat_transfer[0] * self.temperature[:, 0] + + (1 - self.deep_ocean_efficacy) + * self.ocean_heat_transfer[n_box - 1] + * (self.temperature[:, n_box - 2] - self.temperature[:, n_box - 1]) + ) + self.ocean_heat_content_change = np.cumsum( + self.toa_imbalance + * self.timestep + * earth_radius**2 + * 4 + * np.pi + * seconds_per_year + ) + + +def multi_ebm( + configs, + ocean_heat_capacity, + ocean_heat_transfer, + deep_ocean_efficacy, + stochastic_run, + sigma_eta, + sigma_xi, + gamma_autocorrelation, + seed, + use_seed, + forcing_4co2, + timestep, + timebounds, +): + """Create several instances of the EnergyBalanceModel. + + This allows efficient parallel implementation in FaIR. + We have to use a for loop in this function as is does not look like the linear + algebra functions in scipy are naturally parallel. + + Parameters + ---------- + configs : list + A named list of climate configurations. + ocean_heat_capacity : `np.ndarray` + Ocean heat capacity of each layer (top first), W m-2 yr K-1 + ocean_heat_transfer : `np.ndarray` + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : float + efficacy of deepest ocean layer. See e.g. [Geoffroy2013]_. + stochastic_run : bool + Activate the stochastic variability component from [Cummins2020]_. + sigma_eta : float + Standard deviation of stochastic forcing component from [Cummins2020]_. + sigma_xi : float + Standard deviation of stochastic disturbance applied to surface + layer. See [Cummins2020]_. + gamma_autocorrelation : float + Stochastic forcing continuous-time autocorrelation parameter. + See [Cummins2020]_. + seed : int or None + Random seed to use for stochastic variability. + use_seed : bool + Whether or not to use the random seed. + forcing_4co2 : float + effective radiative forcing from a quadrupling of atmospheric + CO2 concentrations above pre-industrial. + timestep : float + Time interval of the model (yr) + timebounds : np.ndarray + Vector representing the time snapshots to calculate temperature on. + """ + n_configs = ocean_heat_capacity.shape[0] + n_layers = ocean_heat_capacity.shape[1] + n_timebounds = len(timebounds) + ebms = xr.Dataset( + { + "eb_matrix_d": ( + ["config", "eb_dim0", "eb_dim1"], + np.ones((n_configs, n_layers + 1, n_layers + 1)) * np.nan, + ), + "forcing_vector_d": ( + ["config", "eb_dim0"], + np.ones((n_configs, n_layers + 1)) * np.nan, + ), + "stochastic_d": ( + ["timebounds", "config", "eb_dim0"], + np.ones((n_timebounds, n_configs, n_layers + 1)) * np.nan, + ), + "ecs": (["config"], np.ones(n_configs) * np.nan), + "tcr": (["config"], np.ones(n_configs) * np.nan), + }, + coords={ + "timebounds": timebounds, + "config": configs, + "eb_dim0": np.arange(-1, n_layers), + "eb_dim1": np.arange(-1, n_layers), + }, + ) + + for i_conf, config in enumerate(configs): + ebm = EnergyBalanceModel( + ocean_heat_capacity=ocean_heat_capacity[i_conf, :], + ocean_heat_transfer=ocean_heat_transfer[i_conf, :], + deep_ocean_efficacy=deep_ocean_efficacy[i_conf], + stochastic_run=stochastic_run[i_conf], + sigma_eta=sigma_eta[i_conf], + sigma_xi=sigma_xi[i_conf], + gamma_autocorrelation=gamma_autocorrelation[i_conf], + seed=seed.data[i_conf] if use_seed[i_conf] else None, + forcing_4co2=forcing_4co2[i_conf], + timestep=timestep, + n_timesteps=n_timebounds, + ) + ebms["eb_matrix_d"].loc[dict(config=config)] = ebm.eb_matrix_d + ebms["forcing_vector_d"].loc[dict(config=config)] = ebm.forcing_vector_d + ebms["stochastic_d"].loc[dict(config=config)] = ebm.stochastic_d + ebm.emergent_parameters() + ebms["ecs"].loc[dict(config=config)] = ebm.ecs + ebms["tcr"].loc[dict(config=config)] = ebm.tcr + + return ebms + + +def step_temperature(state_old, eb_matrix_d, forcing_vector_d, stochastic_d, forcing): + """Advance parallel energy balance models forward one timestep. + + Parameters + ---------- + state_old : np.ndarray + stacked arrays of forcing and temperature of layers in previous timestep + eb_matrix_d : np.ndarray + stacked discretised energy balance matrices + forcing_vector_d : np.ndarray + stacked discretised forcing vectors + _stochastic_d : np.ndarray + stacked matrices of stochastic internal variability + forcing: np.ndarray + stacked vectors of [effective] radiative forcing + + Returns + ------- + state_new : np.ndarray + stacked arrays of forcing and temperature of layers in this timestep + """ + state_new = ( + (eb_matrix_d[0, ...] @ state_old[0, ..., None])[..., 0] + + forcing_vector_d[0, ...] * forcing[0, ..., 0, None] + + stochastic_d[0, ...] + ) + + return state_new + + +def calculate_toa_imbalance_postrun( + state, + forcing, + ocean_heat_transfer, + deep_ocean_efficacy, +): + """Calculate top of atmosphere energy imbalance. + + The calculation is performed after the scenario has been run to avoid + looping, since no dynamic state changes affect the calculation. + + Parameters + ---------- + state : np.ndarray + stacked arrays of forcing and temperature of layers across the run + forcing : np.ndarray + stacked arrays of [effective] radiative forcing across the run + ocean_heat_transfer : np.ndarray + Heat exchange coefficient between ocean layers (top first). The + first element of this array is akin to the climate feedback + parameter, with the convention that stabilising feedbacks are + positive (opposite to most climate sensitivity literature). + W m-2 K-1 + deep_ocean_efficacy : np.ndarray + efficacy of deepest ocean layer. + + Returns + ------- + toa_imbalance : np.ndarray + Top of atmsophere energy imbalance. + """ + toa_imbalance = ( + forcing + - ocean_heat_transfer[..., 0] * state[..., 1] + + (1 - deep_ocean_efficacy) + * ocean_heat_transfer[..., -1] + * (state[..., -2] - state[..., -1]) + ) + return toa_imbalance diff --git a/modules/ebm3/sterodynamics/pipeline.global.yml b/modules/ebm3/sterodynamics/pipeline.global.yml new file mode 100644 index 00000000..4872e770 --- /dev/null +++ b/modules/ebm3/sterodynamics/pipeline.global.yml @@ -0,0 +1,21 @@ +project: + task1: + executable: "python3" + python_dependencies: "numpy scipy netCDF4 pyyaml h5py pandas" + script: "ebm3_thermalexpansion_project.py" + options: + - "scenario" + - "climate_data_file" + - "pipeline_id" + - "nsamps" + - "seed" + - "scenario" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_and_extract_input_data: + - "ebm3_project_data.tgz" + global_total_files: + - "%PIPELINE_ID%_globalsl.nc" diff --git a/modules/ebm3/sterodynamics/pipeline.yml b/modules/ebm3/sterodynamics/pipeline.yml new file mode 100644 index 00000000..cdf00fad --- /dev/null +++ b/modules/ebm3/sterodynamics/pipeline.yml @@ -0,0 +1,57 @@ +project: + task1: + executable: "python3" + python_dependencies: "numpy scipy netCDF4 pyyaml h5py pandas" + script: "ebm3_thermalexpansion_project.py" + options: + - "scenario" + - "climate_data_file" + - "pipeline_id" + - "nsamps" + - "seed" + - "scenario" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + upload_and_extract_input_data: + - "ebm3_project_data.tgz" + global_total_files: + - "%PIPELINE_ID%_globalsl.nc" + +postprocess: + task1: + executable: "python3" + python_dependencies: "scikit-learn tqdm" + script: "ebm3_oceandynamics_project.py" + options: + - "scenario" + - "pipeline_id" + - "nsamps" + - "seed" + - "climate_data_file" + - "baseyear" + - "pyear_start" + - "pyear_end" + - "pyear_step" + - "seed" + - "locationfile" + - "rfmip" + - "params" + - "zosdir" + upload_input_data: + - "energy_balance_model.py" + copy_input_data: + project: + task1: + - "%PIPELINE_ID%_globalsl.nc" + upload_and_extract_input_data: + - "ebm3_project_data.tgz" + - "tlm_sterodynamics_cmip6_data.tgz" + copy_shared_data: + - '$SHARED/location.lst' + download_output_data: + - "%PIPELINE_ID%_quantiles.nc" + local_total_files: + - "%PIPELINE_ID%_localsl.nc" From dfabc79625741d601c9944d63f422be10f3797bb Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 30 Mar 2026 21:11:28 -0700 Subject: [PATCH 08/15] Remove old ebm3/thermalexpansion. Update ebm3/sterodynamics experiment and pipeline for new names and input data --- experiments/ebm3.global/config.yml | 10 +- modules-data/modules-data.urls.txt | 2 + modules/ebm3/sterodynamics/pipeline.yml | 8 +- .../ebm3_thermalexpansion_postprocess.py | 490 ------------------ .../ebm3_thermalexpansion_project.py | 164 ------ .../ebm3/thermalexpansion/pipeline.global.yml | 21 - modules/ebm3/thermalexpansion/pipeline.yml | 54 -- 7 files changed, 9 insertions(+), 740 deletions(-) delete mode 100644 modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py delete mode 100644 modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py delete mode 100644 modules/ebm3/thermalexpansion/pipeline.global.yml delete mode 100644 modules/ebm3/thermalexpansion/pipeline.yml diff --git a/experiments/ebm3.global/config.yml b/experiments/ebm3.global/config.yml index 044f60fa..9670377e 100644 --- a/experiments/ebm3.global/config.yml +++ b/experiments/ebm3.global/config.yml @@ -11,11 +11,7 @@ climate_step: fair2: module_set: "fair2" module: "climate" - pipeline_file: 'pipeline.yml' - # input_data_file: - # - "ebm3.global.fair2.fair2.climate_ohc.nc" - # - "ebm3.global.fair2.fair2.climate_gsat.nc" - # - "ebm3.global.fair2.fair2.climate_climate.nc" + pipeline_file: 'pipeline.global.yml' climate_output_data: - "ebm3.global.fair2.fair2.climate_ohc.nc" - "ebm3.global.fair2.fair2.climate_gsat.nc" @@ -24,6 +20,6 @@ climate_step: generates_climate_output: true sealevel_step: - ebm3: + ocean: module_set: "ebm3" - module: "thermalexpansion" \ No newline at end of file + module: "sterodynamics" \ No newline at end of file diff --git a/modules-data/modules-data.urls.txt b/modules-data/modules-data.urls.txt index 38667d38..9ca15dfd 100644 --- a/modules-data/modules-data.urls.txt +++ b/modules-data/modules-data.urls.txt @@ -30,3 +30,5 @@ https://zenodo.org/records/11506798/files/fair2_climate_project_data.tgz https://zenodo.org/records/11506798/files/ebm3_thermal_expansion_data.tgz https://zenodo.org/records/18199757/files/oelsmann24_vlm_data.tar.gz https://zenodo.org/records/15975768/files/caron18_GIA_preprocess_data.tgz +https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_data.tgz +https://zenodo.org/records/16024082/files/ebm3_project_data.tgz \ No newline at end of file diff --git a/modules/ebm3/sterodynamics/pipeline.yml b/modules/ebm3/sterodynamics/pipeline.yml index cdf00fad..95cb4d8b 100644 --- a/modules/ebm3/sterodynamics/pipeline.yml +++ b/modules/ebm3/sterodynamics/pipeline.yml @@ -16,7 +16,7 @@ project: - "pyear_step" - "seed" upload_and_extract_input_data: - - "ebm3_project_data.tgz" + - "./modules-data/ebm3/ebm3_project_data.tgz" global_total_files: - "%PIPELINE_ID%_globalsl.nc" @@ -47,10 +47,10 @@ postprocess: task1: - "%PIPELINE_ID%_globalsl.nc" upload_and_extract_input_data: - - "ebm3_project_data.tgz" - - "tlm_sterodynamics_cmip6_data.tgz" + - "./modules-data/ebm3/ebm3_project_data.tgz" + - "./modules-data/ebm3/tlm_sterodynamics_cmip6_data.tgz" copy_shared_data: - - '$SHARED/location.lst' + - "$SHARED/location.lst" download_output_data: - "%PIPELINE_ID%_quantiles.nc" local_total_files: diff --git a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py b/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py deleted file mode 100644 index bbb87b05..00000000 --- a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_postprocess.py +++ /dev/null @@ -1,490 +0,0 @@ -import numpy as np -import os -import pandas as pd -import sys -import argparse -from fair import energy_balance_model as ebm3 -import csv -import xarray as xr -from sklearn import linear_model -import re -import time -import netCDF4 -import h5py -import scipy -import pickle -from tqdm import tqdm - -""" -Created on Fri Jun 7 11:24:19 2024 - -For any model you have up tp 2300, you want a set of scaling coefficients -for ssp126 and ssp585, as those thwo a typically the ones available up to 2300. -Those scenarios are used to get a pattern of regression coefficient used for -SSP interpolation and probabilistic projections using FAIR. - -@author: vmalagonsantos -""" - -# functions -def fetch_erfs_from_rcmip(path,scenarios): - ''' - function to read scenario-dependent ERF from RCMIP csv - - Input parameters: - path = path to RCMIP csv - scenarios = desired ssp's - - Output: Adds results of extreme sea-level analysis for all stations to the - station data dictionary. - - erfs = effective radiative forcing per scenario and year - erfyears = corresponding years (1750-2500) - ''' - - erfyears=np.arange(1750,2501) - ssp_idx = {'ssp119':212,'ssp126':231,'ssp245':308,'ssp370':59,'ssp585':404} #location in table - - with open(path) as csv_file: - csv_reader = csv.reader(csv_file) - rows = list(csv_reader) - - erfs = np.empty((len(scenarios),len(erfyears))) - - for s,scen in enumerate(scenarios): - try: - erfs[s,:] = rows[ssp_idx[scen]][7:] - except: - continue - - return erfs, erfyears - -def angd(lat0, lon0, qlat, qlon): - - # Convert the input from degrees to radians - (lat0, lon0) = np.radians((lat0, lon0)) - (qlat, qlon) = np.radians((qlat, qlon)) - - # Calculate the angle between the vectors - temp = np.arctan2(np.sqrt((np.cos(qlat)*np.sin(qlon-lon0))**2 + \ - (np.cos(lat0)*np.sin(qlat) - np.sin(lat0)*np.cos(qlat) * np.cos(qlon-lon0))**2),\ - (np.sin(lat0)*np.sin(qlat) + np.cos(lat0)*np.cos(qlat)*np.cos(qlon-lon0))) - - # Convert the results from radians to degrees and return - return(np.degrees(temp)) - -def NearestPoint(qlat, qlon, lats, lons, tol = None): - - # Get the distance between the query point and all the possible points - dist = angd(lats, lons, qlat, qlon) - - # Which is the closest point - nearest_idx = np.argmin(dist) - - # Is the point within the tolerance? - if isinstance(tol, (int, float)): - if dist[nearest_idx] > tol: - return(None) - - return(nearest_idx) - - -def NearestPoints(qlats, qlons, lats, lons, tol): - - if len(qlats) != len(qlons): - raise Exception("Query lats ({}) and lons ({}) differ in length".format(len(qlats), len(qlons))) - - idx = map(lambda qlat,qlon: NearestPoint(qlat, qlon, lats, lons, tol), qlats, qlons) - - return(list(idx)) - -def ReadLocationFile(location_file): - - # Initialize variables to hold data and site information - names = [] - ids = [] - lats = [] - lons = [] - - # Compile the regex for finding commented lines - comment_regex = re.compile(r'^#') - - # Open the rate file - with open(location_file, 'r') as f: - - # Loop over the lines of the file - for line in f: - - # Skip commented lines - if re.search(comment_regex, line): - continue - - # Split the line into components - (this_name, this_id, this_lat, this_lon) = line.split("\t") - - # Store the information - names.append(this_name) - ids.append(int(this_id)) - lats.append(float(this_lat)) - lons.append(float(this_lon)) - - # Cast everything as numpy arrays - names = np.array(names) - ids = np.array(ids) - lats = np.array(lats) - lons = np.array(lons) - - # Return variables - return(names, ids, lats, lons) - - -# Smooth ZOS and ZOSTOGA over 19 year smoothing window -def nanSmooth(x, w=19): - idx = np.flatnonzero(~np.isnan(x)) - temp = x - if len(idx) > 0: - temp[idx] = Smooth(x[idx], w) - return(temp) - -def Smooth(x, w=19): - out0 = np.convolve(x, np.ones(w,dtype='double'), 'valid')/w - r = np.arange(1,w-1,2, dtype="double") - start = np.cumsum(x[:w-1])[::2]/r - stop = (np.cumsum(x[:-w:-1])[::2]/r)[::-1] - y = np.concatenate((start, out0, stop)) - return(y) - -def emb3_thermalexpansion_postprocess(scenario, pipeline_id, nsamps, seed, pyear_start, pyear_end, pyear_step, locationfile, baseyear, climate_data_file, rfmip, params, zosdir): - targyears = np.arange(baseyear,2301) # for regression - projyears = np.arange(pyear_start,pyear_end+1, pyear_step) - (_, site_ids, site_lats, site_lons) = ReadLocationFile(locationfile) - - # get temperature from FaIR simulations: INPUT FROM CLIMATE STEP. We need both gsat and oceantemp - cf = xr.open_dataset(climate_data_file, group=scenario, engine='netcdf4', ) - gsat = cf['surface_temperature'].sel(years=projyears) - cf['surface_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - otemp = cf['deep_ocean_temperature'].sel(years=projyears) - cf['deep_ocean_temperature'].sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - - # INPUT temperature file from ebm3 global - gte_file = f'{pipeline_id}_globalsl.nc' - - # INPUT get models and parameters - paramdir = params - - #forcing for ebm - scenarios = ['ssp126', 'ssp585'] # NOT INPUT, these two are needed for SSP interpolation in 2300 projections - rfmipfile = rfmip #path # RFMIP FILE IS INPUT, NEW DATA - erfs,erfyears = fetch_erfs_from_rcmip(rfmipfile, scenarios) #get ERF timeseries - - ebm_param = pd.read_csv(paramdir) - - # find available zos simulations - ignoring hidden files - model_has_zos = [f for f in os.listdir(zosdir) if not f.startswith('.')] - - ## FITTING - slopes = [] - intercept = [] - - for m, model in enumerate(model_has_zos): - # model = model_has_zos[0] - - # model = 'MRI-ESM2-0' - - # find matching ebm parameters to model - model_param = ebm_param[ebm_param["model"].str.contains(model)] - - modeldir = zosdir + model - # keep those up to 2300 or 2500 - zos_runs2_2300 = [f for f in os.listdir(modeldir) if '2300' in f or '2500' in f] # regex would be better here - if not zos_runs2_2300: - print('Model ' + model + ' only runs to 2100. Skipping and looking for longer simulations' ) - - continue - - print('Estimating parameters for model: ' + model) - - - # get variants with 2300 simulation - if len(zos_runs2_2300) == 1: - variant_2300 = [f for f in re.split('_', zos_runs2_2300[0]) if f.startswith('r')] - model_var_param = model_param[ebm_param["run"].str.contains(variant_2300[0])] - else: - for v, filename in enumerate(zos_runs2_2300): - variant_2300 = [f for f in re.split('_', zos_runs2_2300[v]) if f.startswith('r')] - # get EBM parameters for both model and variant - model_var_param = model_param[ebm_param["run"].str.contains(variant_2300[0])] - if not model_var_param.empty: - continue - - c1 = np.array(model_var_param.get('C1')) - c2 = np.array(model_var_param.get('C2')) - c3 = np.array(model_var_param.get('C3')) - k1 = np.array(model_var_param.get('kappa1')) - k2 = np.array(model_var_param.get('kappa2')) - k3 = np.array(model_var_param.get('kappa3')) - e = np.array(model_var_param.get('epsilon')) - f4 = np.array(model_var_param.get('F_4xCO2')) - se = np.array(model_var_param.get('sigma_eta')) - sx = np.array(model_var_param.get('sigma_xi')) - gamma = np.array(model_var_param.get('gamma')) - nit = np.array(model_var_param.get('nit')) - - - temp3 = np.empty((len(erfyears), 3, len(scenarios))) - - for s, scen in enumerate(scenarios): - - #Model setup 3lm - ebm_obj = ebm3.EnergyBalanceModel( - ocean_heat_capacity = [c1[0], c2[0], c3[0]], - ocean_heat_transfer =[k1[0], k2[0], k3[0]], ## change this placeholder value - deep_ocean_efficacy = e[0], - forcing_4co2 = f4[0], - stochastic_run=False, - sigma_eta = se[0], - sigma_xi = sx[0], - gamma_autocorrelation = gamma[0], - seed = None, - timestep = 1, - n_timesteps = 2500-1750, - ) - - # Run model for different scenarios - ebm_obj.add_forcing( erfs[s,:], erfyears) - ebm_obj.run() - temp3[:,:,s] = ebm_obj.temperature - - attrs={'description': 'surface and deep temperature obtained by EBM-3LM', - 'contact':'Victor Malagon Santos, victor.malagon.santos@nioz.nl'} - - temp3_xr = xr.Dataset({'temperature':(['years','layers', 'scenarios'], temp3)}, - coords={'layers':['surface', 'inter', 'deep'], 'years': erfyears, 'scenarios': scenarios}, attrs=attrs) - - # reference to base year - temp3_xr = temp3_xr - temp3_xr.sel(years=baseyear) - # .sel(years=np.arange(baseyear-9,baseyear+10)).mean(dim='years') - - # get zos, historical + ssp - model_dir = zosdir + model - - hist_file = [f for f in os.listdir(model_dir) if 'hist' in f and variant_2300[0] in f] - ssp1_file = [f for f in os.listdir(model_dir) if 'ssp1' in f and variant_2300[0] in f and 'ssp126' in f] - ssp5_file = [f for f in os.listdir(model_dir) if 'ssp5' in f and variant_2300[0] in f and 'ssp585' in f] - - if not ssp1_file or not ssp5_file: - continue - - hist_xr = xr.open_dataset(model_dir + '/' + hist_file[0],decode_times=False) - ssp1_xr = xr.open_dataset(model_dir + '/' + ssp1_file[0],decode_times=False) - ssp5_xr = xr.open_dataset(model_dir + '/' + ssp5_file[0],decode_times=False) - - if ssp1_xr['time'][0].values != 60280.5: - ssp1_xr.coords['time'] = ((ssp1_xr['time'] + 60280.5).astype('float')) - ssp5_xr.coords['time'] = ((ssp5_xr['time'] + 60280.5).astype('float')) - - # cocatenate - zos1 = xr.concat([hist_xr, ssp1_xr], dim='time') - zos5 = xr.concat([hist_xr, ssp5_xr], dim='time') - - #% create new coordinates for year - zos1.coords['years'] = ((zos1['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR - zos5.coords['years'] = ((zos5['time'] // 365.25)+1850).astype('int') # REMOVE HARDCODED STARTING YEAR - - # Calculate annual means - zos1_am = zos1.groupby('years').mean(dim='time') - zos1_am = zos1_am.where(zos1_am['zos'] < 99999) # replace land values with nans - zos5_am = zos5.groupby('years').mean(dim='time') - zos5_am = zos5_am.where(zos5_am['zos'] < 99999) # replace land values with nans - - # reference zos - # baseyear_idx = np.flatnonzero(datayears == baseyear) - # S = np.apply_along_axis(lambda z, idx: z - z[idx], axis=0, arr=sZOS, idx=baseyear_idx) - zos1_am = zos1_am.sel(years=targyears) - zos1_am.sel(years=baseyear) - zos5_am = zos5_am.sel(years=targyears) - zos5_am.sel(years=baseyear) - - #% fitting stage - nlon = len(zos1_am['lon']) - nlat = len(zos1_am['lat']) - - - Ts = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='surface').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='surface').drop_vars('scenarios')),dim='years') - - Ti = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='inter').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='inter').drop_vars('scenarios')),dim='years') - - Td = xr.concat((temp3_xr.sel(years=targyears).sel(scenarios=['ssp126']).sel(layers='deep').drop_vars('scenarios'), - temp3_xr.sel(years=targyears).sel(scenarios=['ssp585']).sel(layers='deep').drop_vars('scenarios')),dim='years') - - - zos1 = np.array(zos1_am['zos'].sel(years=targyears)) - zos5 = np.array(zos5_am['zos'].sel(years=targyears)) - - slope = np.zeros((3, nlat, nlon)) - slope[:] = np.nan - smoothwin = 19; - - reg = linear_model.LinearRegression() - # for i in tqdm(range(nlat)): - # sleep(3) - for i in tqdm(range(nlat)): - for j in range(nlon): - - y = np.concatenate((Smooth(zos1[:,i,j],w=smoothwin), Smooth(zos5[:,i,j],w=smoothwin))) # reducing varability - x1 = np.array(Ts['temperature'][:]).flatten() - x2 = np.array(Ti['temperature'][:]).flatten() - x3 = np.array(Td['temperature'][:]).flatten() - x = np.transpose(np.stack((x1,x2,x3))) - - try: - - reg.fit(x, y) - slope[0,i,j] = reg.coef_[0] # slope coefficient of T surface - slope[1,i,j] = reg.coef_[1] # slope coefficient of T inter - slope[2,i,j] = reg.coef_[2] # slope coefficient of T deep - # intercept[m,i,j] = reg.intercept_ - - except ValueError: - slope[0,i,j] = np.nan - slope[1,i,j] = np.nan - slope[2,i,j] = np.nan - # intercept[m,i,j] = np.nan - - slopes.append(slope) # MAYBE SAVE TO A PICKLE, TO BE PROVIDED IN THE NEXT TASK? - - lons, lats = np.meshgrid(zos1_am['lon'], zos1_am['lat']) - lats = lats.flatten() - lons = lons.flatten() - - site_ids_map = np.array(NearestPoints(site_lats, site_lons, lats, lons, tol=None)) - - #nsamps = 100 - samples = np.array(gsat['samples']) - - # RESAMPLE SLOPE PARAMETERS - nsims=len(slopes) - - rng = np.random.default_rng(seed) - if nsamps > nsims: - run_idx = np.arange(nsims) - sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) - else: - run_idx = rng.choice(nsims, nsamps, nsamps>nsims) - sample_idx = np.arange(nsamps) - - - ## need to figure out of to do this resamples properly, so from now on this is dummy data - slopes_resampled = [slopes[i] for i in sample_idx] - # aa = np.dstack(slopes_resampled) - - dsl = np.empty((nsamps, len(projyears), len(site_ids))) - dsl[:] = np.nan - - for sample in tqdm(samples): - slope = slopes_resampled[sample] - print(sample) - - # prepare slopes - slope_s = slope[0,:,:].flatten() - slope_s = np.array([slope_s[x] for x in site_ids_map]) - - slope_i = slope[1,:,:].flatten() - slope_i = np.array([slope_i[x] for x in site_ids_map]) - - slope_d = slope[2,:,:].flatten() - slope_d = np.array([slope_d[x] for x in site_ids_map]) - - # obtain fair temepratures - Tfs = np.array(gsat.sel(samples=sample)) # surface temperature from fair - Tfi = np.array(otemp.sel(layers=1).sel(samples=sample)) # intermediate temperature from fair - Tfd = np.array(otemp.sel(layers=2).sel(samples=sample)) # deep temperature from fair - - - #% project - dsl[sample,:,:] = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) - - # dsl = np.multiply.outer(Tfs, slope_s) + np.multiply.outer(Tfi, slope_i) + np.multiply.outer(Tfd, slope_d) - - ncvar_attributes = {"description": "Dynamic Sea Level"} - - nc_missing_value = np.nan - # Generate the output xarray - dsl_xr = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), dsl, {"units":"mm", "missing_value":nc_missing_value}), - "lat": (("locations"), site_lats), - "lon": (("locations"), site_lons)}, - coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) - - - # make sure projections have 0 mean. This weighted average only works for regular grids. - # Must be edited to accomodate other grids if zos is not in 1x1. - - weights = np.cos(np.deg2rad(dsl_xr.lat)) - weights.name = "weights" - - dsl_xr_weighted = dsl_xr.weighted(weights) - dsl_xr_weighted - weighted_mean = dsl_xr_weighted.mean(("locations")) - - - # Appply 0mean correction to dsl, and add GTE to get sterodynamics - - gte = xr.open_dataset(gte_file) - - - gte = gte['sea_level_change'].values[:,:,0] - sdsl = dsl*1000 - np.repeat(np.array(weighted_mean['sea_level_change'])[:, :, np.newaxis], len(site_ids), axis=2) + np.repeat(gte[0:nsamps, :, np.newaxis], len(site_ids), axis=2) - - ncvar_attributes = {"description": "Local SLR contributions from thermal expansion and dynamic sea-level using EBM3", - "history": "Created " + time.ctime(time.time()), - "source": "SLR Framework: PROTECT 2300", - "scenario": scenario, # change to scenario later - "baseyear": baseyear} - - nc_missing_value = np.nan - # Generate the output xarray - local_out = xr.Dataset({"sea_level_change": (("samples", "years", "locations"), sdsl, {"units":"mm", "missing_value":nc_missing_value}), - "lat": (("locations"), site_lats), - "lon": (("locations"), site_lons)}, - coords={"years": projyears, "locations": site_ids, "samples": np.arange(nsamps)}, attrs=ncvar_attributes) - # Write these samples to a temporary netcdf file - local_out.to_netcdf("{0}_localsl.nc".format(pipeline_id), encoding={"sea_level_change": {"dtype": "f4", "zlib": True, "complevel":4, "_FillValue": nc_missing_value}}) - -if __name__ == '__main__': - - # Initialize the command-line argument parser - parser = argparse.ArgumentParser(description="Run the pre-processing stage for the TLM ocean dynamics workflow",\ - epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") - - # Define the command line arguments to be expected - parser.add_argument('--scenario', help="SSP scenario (i.e ssp585) or temperature target (i.e. tlim2.0win0.25)", default='ssp585') - parser.add_argument('--nsamps', help="Number of samples to generate [default=20000]", default=20000, type=int) - parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) - parser.add_argument('--pyear_start', help="Year for which projections start [default=2000]", default=2000, type=int) - parser.add_argument('--pyear_end', help="Year for which projections end [default=2300]", default=2300, type=int) - parser.add_argument('--pyear_step', help="Step size in years between pyear_start and pyear_end at which projections are produced [default=5]", default=5, type=int) - parser.add_argument('--locationfile', help="File that contains name, id, lat, and lon of points for localization", default="location.lst") - parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) - parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") - parser.add_argument('--climate_data_file',type=str) - parser.add_argument('--rfmip', help='rfmip file',default='rfmip-radiative-forcing-annual-means-v4-0-0.csv') - parser.add_argument('--params', help='CMIP6 Params cvs', default='4xCO2_cummins_ebm3_cmip6.csv') - parser.add_argument('--zosdir',help='Path to CMIP6 ZOS directory', default='cmip6/zos/') - - # Parse the arguments - args = parser.parse_args() - - emb3_thermalexpansion_postprocess(args.scenario, - args.pipeline_id, - args.nsamps, - args.seed, - args.pyear_start, - args.pyear_end, - args.pyear_step, - args.locationfile, - args.baseyear, - args.climate_data_file, - args.rfmip, - args.params, - args.zosdir) - - # Done - sys.exit() \ No newline at end of file diff --git a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py b/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py deleted file mode 100644 index 6b882dbd..00000000 --- a/modules/ebm3/thermalexpansion/ebm3_thermalexpansion_project.py +++ /dev/null @@ -1,164 +0,0 @@ -# GMSLR projection program used for IPCC WG1 AR5 -# Translated from IDL to Python 2.7 by Jonathan Gregory 23.10.19 - -import os -import numpy as np -import netCDF4 as nc -import pandas as pd -import time -from netCDF4 import Dataset -import pickle -import sys -import re -import argparse - -class ProjectionError(Exception): - pass - - -def ebm3_thermalexpansion_project(scenario, climate_data_file, coef_file, params_file, pyear_start, pyear_end, pyear_step, nsamps, pipeline_id, seed, baseyear): - # constants - a = 6.37*1e6 - earth_area = 4*np.pi*a**2 - - path = os.path.dirname(__file__) - - # heat capacity of each layer from fair2. Needed to compute OHC - fparam = pd.read_csv(params_file) - c1 = fparam['clim_c1'] - c2 = fparam['clim_c2'] - c3 = fparam['clim_c3'] - - nsims=1001 - - rng = np.random.default_rng(seed) - if nsamps > nsims: - run_idx = np.arange(nsims) - sample_idx = rng.choice(nsims, nsamps, nsamps>nsims) - else: - run_idx = rng.choice(nsims, nsamps, nsamps>nsims) - sample_idx = np.arange(nsamps) - - c1 = c1[sample_idx] - c2 = c2[sample_idx] - c3 = c3[sample_idx] - - # create target years array - targyears = np.arange(pyear_start, pyear_end+1, pyear_step) - - #temperature output from fair2. Needed t compute OHC - ds = Dataset(climate_data_file) - gsat = ds[scenario]['surface_temperature'][:] - deepoceant = ds[scenario]['deep_ocean_temperature'][:] - years = ds[scenario]['years'][:] - - # Expansion coefficients. Needed to to convert OCH to global thermal expansion (GTE) - ds = Dataset(coef_file) - include_models = ds['model'][:] - eeh3 = ds['expcoefs'][:] - - # estimate OHC - up = gsat*np.array(c1) - mid = deepoceant[:,:,0]*np.array(c2) - low = deepoceant[:,:,1]*np.array(c3) - n = np.zeros((len(years), nsamps)) - n[1:,:] = up[1:,:]-up[:-1,:] + mid[1:,:]-mid[:-1,:] + low[1:,:]-low[:-1,:] - ohc_samps = np.cumsum(n*earth_area, 0) * 365*24*3.6*1e3*1e-24 - - # Generate samples assuming normal distribution - rng = np.random.default_rng(seed) - expcoef_samps = rng.normal(loc=np.mean(eeh3), scale=np.std(eeh3), size=(nsamps,1)) - - # Produce the projection samples - gte_samps = ohc_samps * expcoef_samps.flatten() - - # Center these samples on the baseyear - baseyear_idx = np.flatnonzero(years == baseyear) - gte_samps = gte_samps - gte_samps[baseyear_idx,:] - - # Subset the samples for the projection years - targyear_idx = np.isin(years, targyears) - gte_samps = gte_samps[targyear_idx,:] - - # Invert the dimensions of the variable and convert from m to mm - gte_samps *= 1000. - gte_samps = gte_samps.T - - # Save the projections to a pickle - output = {"thermsamps": gte_samps, "targyears": targyears, "baseyear": baseyear, \ - "include_models": include_models, "scenario": scenario} - outfile = open(os.path.join(os.path.dirname(__file__), "{}_projections.pkl".format(pipeline_id)), 'wb') - pickle.dump(output, outfile) - outfile.close() - - # Write the total global projections to a netcdf file - nc_filename = os.path.join(os.path.dirname(__file__), "{0}_globalsl.nc".format(pipeline_id)) - rootgrp = Dataset(nc_filename, "w", format="NETCDF4") - - # Define Dimensions - nyr = len(targyears) - year_dim = rootgrp.createDimension("years", nyr) - samp_dim = rootgrp.createDimension("samples", nsamps) - loc_dim = rootgrp.createDimension("locations", 1) - - # Populate dimension variables - year_var = rootgrp.createVariable("years", "i4", ("years",)) - samp_var = rootgrp.createVariable("samples", "i8", ("samples",)) - loc_var = rootgrp.createVariable("locations", "i8", ("locations",)) - lat_var = rootgrp.createVariable("lat", "f4", ("locations",)) - lon_var = rootgrp.createVariable("lon", "f4", ("locations",)) - - # Create a data variable - samps = rootgrp.createVariable("sea_level_change", "f4", ("samples", "years", "locations"), zlib=True, complevel=4) - - # Assign attributes - rootgrp.description = "Global SLR contribution from Thermal Expansion according to Three-Layer Model workflow" - rootgrp.history = "Created " + time.ctime(time.time()) - rootgrp.source = "FACTS: {0}".format(pipeline_id) - rootgrp.scenario = scenario - rootgrp.baseyear = baseyear - rootgrp.comment = "Included Models: " + ",".join([str(x) for x in include_models]) - samps.units = "mm" - - # Put the data into the netcdf variables - year_var[:] = targyears - samp_var[:] = np.arange(nsamps) - samps[:,:,:] = gte_samps[:,:,np.newaxis] - lat_var[:] = np.inf - lon_var[:] = np.inf - loc_var[:] = -1 - - # Close the netcdf - rootgrp.close() - - return(0) - - - -if __name__ == '__main__': - - # Initialize the command-line argument parser - parser = argparse.ArgumentParser(description="Run the glacier projection stage for the AR5 SLR projection workflow",\ - epilog="Note: This is meant to be run as part of the Framework for the Assessment of Changes To Sea-level (FACTS)") - - # Define the command line arguments to be expected - parser.add_argument('--scenario', help="SSP scenario (i.e ssp585) or temperature target (i.e. tlim2.0win0.25)", default='ssp585') - parser.add_argument('--climate_data_file', help="NetCDF4/HDF5 file containing surface temperature data", type=str) - parser.add_argument('--params_file', help='Full path to calibrated constraints params file', default='calibrated_constrained_parameters.csv') - parser.add_argument('--coef_file', help='Full path to expansion coefficient file', default='scmpy3LM_RCMIP_CMIP6calpm_n18_expcoefs.nc') - parser.add_argument('--nsamps', help="Number of samples to generate [default=1000]", default=1000, type=int) - parser.add_argument('--pyear_start', help="Projection year start [default=2020]", default=2020, type=int) - parser.add_argument('--pyear_end', help="Projection year end [default=2100]", default=2150, type=int) - parser.add_argument('--pyear_step', help="Projection year step [default=10]", default=10, type=int) - parser.add_argument('--baseyear', help="Base year to which slr projections are centered", type=int, default=2005) - parser.add_argument('--seed', help="Seed value for random number generator [default=1234]", default=1234, type=int) - parser.add_argument('--pipeline_id', help="Unique identifier for this instance of the module") - - - # Parse the arguments - args = parser.parse_args() - - # Run the projection process on the files specified from the command line argument - ebm3_thermalexpansion_project(args.scenario, args.climate_data_file, args.coef_file, args.params_file, args.pyear_start, args.pyear_end, args.pyear_step, args.nsamps, args.pipeline_id, args.seed, args.baseyear) - - exit() diff --git a/modules/ebm3/thermalexpansion/pipeline.global.yml b/modules/ebm3/thermalexpansion/pipeline.global.yml deleted file mode 100644 index 65518ae1..00000000 --- a/modules/ebm3/thermalexpansion/pipeline.global.yml +++ /dev/null @@ -1,21 +0,0 @@ -project: - task1: - executable: "python3" - python_dependencies: "numpy scipy netCDF4 pyyaml h5py pandas" - script: "ebm3_thermalexpansion_project.py" - options: - - "scenario" - - "climate_data_file" - - "pipeline_id" - - "nsamps" - - "seed" - - "scenario" - - "baseyear" - - "pyear_start" - - "pyear_end" - - "pyear_step" - - "seed" - upload_and_extract_input_data: - - "ebm3_thermal_expansion_data.tgz" - global_total_files: - - "%PIPELINE_ID%_globalsl.nc" diff --git a/modules/ebm3/thermalexpansion/pipeline.yml b/modules/ebm3/thermalexpansion/pipeline.yml deleted file mode 100644 index bb97291b..00000000 --- a/modules/ebm3/thermalexpansion/pipeline.yml +++ /dev/null @@ -1,54 +0,0 @@ -project: - task1: - executable: "python3" - python_dependencies: "numpy scipy netCDF4 pyyaml h5py pandas" - script: "ebm3_thermalexpansion_project.py" - options: - - "scenario" - - "climate_data_file" - - "pipeline_id" - - "nsamps" - - "seed" - - "scenario" - - "baseyear" - - "pyear_start" - - "pyear_end" - - "pyear_step" - - "seed" - upload_and_extract_input_data: - - "ebm3_thermal_expansion_data.tgz" - global_total_files: - - "%PIPELINE_ID%_globalsl.nc" - -postprocess: - task1: - executable: "python3" - python_dependencies: "scikit-learn tqdm" - script: "ebm3_thermalexpansion_postprocess.py" - options: - - "scenario" - - "pipeline_id" - - "nsamps" - - "seed" - - "climate_data_file" - - "baseyear" - - "pyear_start" - - "pyear_end" - - "pyear_step" - - "seed" - - "locationfile" - - "rfmip" - - "params" - - "zosdir" - copy_input_data: - project: - task1: - - "%PIPELINE_ID%_globalsl.nc" - upload_and_extract_input_data: - - "ebm3_thermal_expansion_postprocessing_data.tgz" - - "tlm_sterodynamics_cmip6_data.tgz" - copy_shared_data: - - '$SHARED/location.lst' - local_total_files: - - "%PIPELINE_ID%_localsl.nc" - From a7c1547934d2432dec97d0e81da54363fdbdf737 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 30 Mar 2026 21:13:40 -0700 Subject: [PATCH 09/15] Update quickstart guide to include warning about tlm and ebm3 same name input data file --- docs/source/quickstart.rst | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst index 76f9a931..34393f9c 100644 --- a/docs/source/quickstart.rst +++ b/docs/source/quickstart.rst @@ -29,6 +29,12 @@ Installing and Using FACTS on a GNU/Linux Workstation wget -P facts/modules-data -i facts/modules-data/modules-data.global_only.urls.txt +.. warning:: + One of the input data files for tlm/sterodynamics (https://zenodo.org/record/7478192/files/tlm_sterodynamics_cmip6_data.tgz) has + the same name as the more recent ebm3/sterodynamics module (https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_data.tgz). + They are not interchangeable, as of April 2026. + Workaround: When you download the ebm3 input data files, put them inside a `./modules-data/ebm3/` directory. The ebm3 pipeline.yml is configured with that path. + 3. Create and activate a Python virtual environment, and install FACTS's Python dependences in it. You can use `venv`, `conda` or `virtualenv` to create your Python virtual environment. See `these instructions `_ for further details. Using `venv`:: python3 -m venv ve3 From 4e447032476cb4336e01c0731697e9a0984994fa Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Mon, 30 Mar 2026 22:21:06 -0700 Subject: [PATCH 10/15] add an ebm3.local experiment example --- experiments/ebm3.local/config.yml | 25 +++++++++++++++++++++++++ 1 file changed, 25 insertions(+) create mode 100644 experiments/ebm3.local/config.yml diff --git a/experiments/ebm3.local/config.yml b/experiments/ebm3.local/config.yml new file mode 100644 index 00000000..23d68c88 --- /dev/null +++ b/experiments/ebm3.local/config.yml @@ -0,0 +1,25 @@ +global-options: + nsamps: 50 + scenario: ssp585 + pyear_start: 2020 + pyear_end: 2150 + pyear_step: 10 + baseyear: 2005 + pipeline_file: "pipeline.yml" + +climate_step: + fair2: + module_set: "fair2" + module: "climate" + pipeline_file: 'pipeline.yml' + climate_output_data: + - "ebm3.local.fair2.fair2.climate_ohc.nc" + - "ebm3.local.fair2.fair2.climate_gsat.nc" + - "ebm3.local.fair2.fair2.climate_climate.nc" + - "ebm3.local.fair2.fair2.climate_oceantemp.nc" + generates_climate_output: true + +sealevel_step: + ocean: + module_set: "ebm3" + module: "sterodynamics" \ No newline at end of file From 753f8a442ac12b2dd68d6ff2ad01dcf785554dda Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Tue, 31 Mar 2026 15:51:44 -0700 Subject: [PATCH 11/15] update urls and ebm3 pipeline for new (non-conflicting) ebm3 input filename --- docs/source/quickstart.rst | 6 ------ modules-data/modules-data.global_only.urls.txt | 2 +- modules-data/modules-data.urls.txt | 2 +- modules/ebm3/sterodynamics/pipeline.yml | 6 +++--- 4 files changed, 5 insertions(+), 11 deletions(-) diff --git a/docs/source/quickstart.rst b/docs/source/quickstart.rst index 34393f9c..76f9a931 100644 --- a/docs/source/quickstart.rst +++ b/docs/source/quickstart.rst @@ -29,12 +29,6 @@ Installing and Using FACTS on a GNU/Linux Workstation wget -P facts/modules-data -i facts/modules-data/modules-data.global_only.urls.txt -.. warning:: - One of the input data files for tlm/sterodynamics (https://zenodo.org/record/7478192/files/tlm_sterodynamics_cmip6_data.tgz) has - the same name as the more recent ebm3/sterodynamics module (https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_data.tgz). - They are not interchangeable, as of April 2026. - Workaround: When you download the ebm3 input data files, put them inside a `./modules-data/ebm3/` directory. The ebm3 pipeline.yml is configured with that path. - 3. Create and activate a Python virtual environment, and install FACTS's Python dependences in it. You can use `venv`, `conda` or `virtualenv` to create your Python virtual environment. See `these instructions `_ for further details. Using `venv`:: python3 -m venv ve3 diff --git a/modules-data/modules-data.global_only.urls.txt b/modules-data/modules-data.global_only.urls.txt index 48956186..58ee43b6 100644 --- a/modules-data/modules-data.global_only.urls.txt +++ b/modules-data/modules-data.global_only.urls.txt @@ -10,4 +10,4 @@ https://zenodo.org/record/7478192/files/larmip_icesheet_project_data.tgz https://zenodo.org/record/7478192/files/ssp_landwaterstorage_preprocess_data.tgz https://zenodo.org/record/7478192/files/tlm_sterodynamics_preprocess_data.tgz https://zenodo.org/records/11506798/files/fair2_climate_project_data.tgz -https://zenodo.org/records/11506798/files/ebm3_thermal_expansion_data.tgz \ No newline at end of file +https://zenodo.org/records/16024082/files/ebm3_project_data.tgz \ No newline at end of file diff --git a/modules-data/modules-data.urls.txt b/modules-data/modules-data.urls.txt index 9ca15dfd..ea2ecb56 100644 --- a/modules-data/modules-data.urls.txt +++ b/modules-data/modules-data.urls.txt @@ -30,5 +30,5 @@ https://zenodo.org/records/11506798/files/fair2_climate_project_data.tgz https://zenodo.org/records/11506798/files/ebm3_thermal_expansion_data.tgz https://zenodo.org/records/18199757/files/oelsmann24_vlm_data.tar.gz https://zenodo.org/records/15975768/files/caron18_GIA_preprocess_data.tgz -https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_data.tgz +https://zenodo.org/records/16024082/files/tlm_sterodynamics_cmip6_ebm3_data.tgz https://zenodo.org/records/16024082/files/ebm3_project_data.tgz \ No newline at end of file diff --git a/modules/ebm3/sterodynamics/pipeline.yml b/modules/ebm3/sterodynamics/pipeline.yml index 95cb4d8b..06586c6d 100644 --- a/modules/ebm3/sterodynamics/pipeline.yml +++ b/modules/ebm3/sterodynamics/pipeline.yml @@ -16,7 +16,7 @@ project: - "pyear_step" - "seed" upload_and_extract_input_data: - - "./modules-data/ebm3/ebm3_project_data.tgz" + - "ebm3_project_data.tgz" global_total_files: - "%PIPELINE_ID%_globalsl.nc" @@ -47,8 +47,8 @@ postprocess: task1: - "%PIPELINE_ID%_globalsl.nc" upload_and_extract_input_data: - - "./modules-data/ebm3/ebm3_project_data.tgz" - - "./modules-data/ebm3/tlm_sterodynamics_cmip6_data.tgz" + - "ebm3_project_data.tgz" + - "tlm_sterodynamics_cmip6_ebm3_data.tgz" copy_shared_data: - "$SHARED/location.lst" download_output_data: From 21d8591f2069424c22de7c69da4c8c9d9860d193 Mon Sep 17 00:00:00 2001 From: pkjr Date: Mon, 11 May 2026 14:57:20 +0000 Subject: [PATCH 12/15] Fix emulandice dependency bundling failure handling --- modules/emulandice2/emulandice_config.sh | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/modules/emulandice2/emulandice_config.sh b/modules/emulandice2/emulandice_config.sh index 291309a5..0344473b 100755 --- a/modules/emulandice2/emulandice_config.sh +++ b/modules/emulandice2/emulandice_config.sh @@ -1,7 +1,8 @@ #!/bin/bash +set -e -ROOTDIR=`dirname $0` -MODE=$1 +ROOTDIR="$(dirname "$0")" +MODE=${1:-} case $MODE in "clean") From 53376f3aecb2dcb00d9aee881ff716f8a5b6b6fc Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Fri, 15 May 2026 00:23:50 -0700 Subject: [PATCH 13/15] add dask[complete] to emulandice2 pipelines and update emulandice2.protect config to local --- experiments/emulandice2.protect.ssp585/config.yml | 2 +- modules/emulandice2/GrIS/pipeline.global.2300.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/experiments/emulandice2.protect.ssp585/config.yml b/experiments/emulandice2.protect.ssp585/config.yml index 267451e5..1f1c110f 100644 --- a/experiments/emulandice2.protect.ssp585/config.yml +++ b/experiments/emulandice2.protect.ssp585/config.yml @@ -5,7 +5,7 @@ global-options: pyear_end: 2300 #2150 pyear_step: 10 baseyear: 2005 - pipeline_file: "pipeline.global.2300.yml" # "pipeline.yml" + pipeline_file: "pipeline.2300.yml" # "pipeline.yml" climate_step: temperature: diff --git a/modules/emulandice2/GrIS/pipeline.global.2300.yml b/modules/emulandice2/GrIS/pipeline.global.2300.yml index 8e981c11..3d265a37 100644 --- a/modules/emulandice2/GrIS/pipeline.global.2300.yml +++ b/modules/emulandice2/GrIS/pipeline.global.2300.yml @@ -1,7 +1,7 @@ project: task1: - python_dependencies: "xarray dask numpy scipy" + python_dependencies: "xarray dask[complete] numpy scipy" executable: "python3" script_noupload: "emulandice_project.py" arguments: From 9f3ff384c347bf0e558216e81e7cf4754a6adf38 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Fri, 15 May 2026 00:25:45 -0700 Subject: [PATCH 14/15] add dask[complete] to emulandice2 pipelines --- modules/emulandice2/AIS/pipeline.2300.yml | 2 +- modules/emulandice2/GrIS/pipeline.2300.yml | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/emulandice2/AIS/pipeline.2300.yml b/modules/emulandice2/AIS/pipeline.2300.yml index 42c507d8..06ea1919 100644 --- a/modules/emulandice2/AIS/pipeline.2300.yml +++ b/modules/emulandice2/AIS/pipeline.2300.yml @@ -1,7 +1,7 @@ project: task1: - python_dependencies: "xarray dask numpy scipy" + python_dependencies: "xarray dask[complete] numpy scipy" executable: "python3" script_noupload: "emulandice_project.py" arguments: diff --git a/modules/emulandice2/GrIS/pipeline.2300.yml b/modules/emulandice2/GrIS/pipeline.2300.yml index 246bb416..6a332691 100644 --- a/modules/emulandice2/GrIS/pipeline.2300.yml +++ b/modules/emulandice2/GrIS/pipeline.2300.yml @@ -1,7 +1,7 @@ project: task1: - python_dependencies: "xarray dask numpy scipy" + python_dependencies: "xarray dask[complete] numpy scipy" executable: "python3" script_noupload: "emulandice_project.py" arguments: From 41c2016d7206e7f49873586811b83d1661e95188 Mon Sep 17 00:00:00 2001 From: Kelly McCusker Date: Thu, 21 May 2026 13:56:02 -0700 Subject: [PATCH 15/15] Fix bug assigning values to numpy array in oelsmann24 --- modules/oelsmann24/verticallandmotion/vlm_preprocess.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/modules/oelsmann24/verticallandmotion/vlm_preprocess.py b/modules/oelsmann24/verticallandmotion/vlm_preprocess.py index b9598b4f..4014a76e 100644 --- a/modules/oelsmann24/verticallandmotion/vlm_preprocess.py +++ b/modules/oelsmann24/verticallandmotion/vlm_preprocess.py @@ -596,7 +596,7 @@ def interpolate_and_compute_weights(VLM_REC_MERGED,VLM_REC_MERGED_GRW,GIA_VLM,di lat_grid = np.arange((180*4))/4-90 distance_max = 150 - lon_ext = dist2coast['lon'].values + lon_ext = dist2coast['lon'].values.copy() lon_ext[lon_ext>180]=lon_ext[lon_ext>180]-360 dist2coast.assign_coords(lon=dist2coast.lon*0+lon_ext) @@ -631,7 +631,7 @@ def interpolate_and_compute_weights(VLM_REC_MERGED,VLM_REC_MERGED_GRW,GIA_VLM,di w_d = w_d.fillna(0) - lon_prt =GIA_VLM['lon'].values + lon_prt =GIA_VLM['lon'].values.copy() lon_prt[lon_prt>180] = lon_prt[lon_prt>180]-360 GIA_VLM['lon'] = lon_prt