From 310835164ac0c420a54fcb11ced1455f4863ef7f Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 27 Apr 2026 18:37:42 +0200 Subject: [PATCH 1/3] Add CLM history file output around DA update` Calls `hist_init_da_tape` after CLM initialisation to register a DA history tape. Introduces `hist_write_pdaf.F90` with two C-callable entry points (clm_hist_write_da_before/after) that capture the forecast and analysis states as standard CLM history files. `enkfpf.par` input variable `CLM:print_da_hist_file` Co-Authored-By: Claude Sonnet 4.6 --- .../running_tsmp_pdaf/input_enkfpf.md | 12 +++ interface/model/Makefile | 1 + interface/model/clm3_5/enkf_clm_mod.F90 | 1 + interface/model/common/enkf.h | 3 + interface/model/common/read_enkfpar.c | 1 + interface/model/eclm/enkf_clm_5.F90 | 2 + interface/model/eclm/enkf_clm_mod_5.F90 | 1 + interface/model/eclm/hist_write_pdaf.F90 | 86 +++++++++++++++++++ interface/model/wrapper_tsmp.c | 6 ++ 9 files changed, 113 insertions(+) create mode 100644 interface/model/eclm/hist_write_pdaf.F90 diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index 790b4aa09..6032ef959 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md @@ -504,6 +504,17 @@ and include the specifier `update` in the file name. `CLM:print_et`: (integer) Invoke function `write_clm_statistics`. For further information, see source code. Default: `0`. +### CLM:print_da_hist_file ### + +`CLM:print_da_hist_file`: (integer) If set to `1`, CLM history file +snapshots of all fields registered on tape 1 (h0) are written before +and after each DA update using the standard CLM history file +infrastructure. One file per DA cycle is produced for each phase, +written to the run directory and named +`caseid.clm2.da_bef.YYYY-MM-DD-SSSSS.nc` and +`caseid.clm2.da_aft.YYYY-MM-DD-SSSSS.nc`, where the +timestamp is the model time of the DA update. Default: `0`. + ### CLM:statevec_allcol ### `CLM:statevec_allcol`: (integer) Switch for using all SWC columns of a @@ -894,6 +905,7 @@ Default: 0, output turned off. | | `update_swc` | 1 | | | `print_swc` | 0 | | | `print_et` | 0 | +| | `print_da_hist_file` | 0 | | | `statevec_allcol` | 0 | | | `statevec_colmean` | 0 | | | `statevec_only_active` | 0 | diff --git a/interface/model/Makefile b/interface/model/Makefile index 185b906aa..35859efad 100644 --- a/interface/model/Makefile +++ b/interface/model/Makefile @@ -26,6 +26,7 @@ OBJCLM = enkf_clm_mod.o\ OBJCLM5 = enkf_clm_mod_5.o\ mod_clm_statistics_5.o\ print_update_clm_5.o\ + hist_write_pdaf.o\ enkf_clm_5.o\ ## parflow object files diff --git a/interface/model/clm3_5/enkf_clm_mod.F90 b/interface/model/clm3_5/enkf_clm_mod.F90 index 7462a6f1e..615065541 100755 --- a/interface/model/clm3_5/enkf_clm_mod.F90 +++ b/interface/model/clm3_5/enkf_clm_mod.F90 @@ -68,6 +68,7 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc + integer(c_int),bind(C,name="clmprint_da_hist_file") :: clmprint_da_hist_file #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index f943ad903..4b77cd840 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -43,6 +43,8 @@ extern void clm_advance(int *ntstep, int *tstartcycle, int *mype); extern void update_clm(int *tstartcycle, int *mype); #if defined CLMSA extern void print_update_clm(int *ts, int *ttot); +extern void clm_hist_write_da_before(); +extern void clm_hist_write_da_after(); #endif extern void write_clm_statistics(int *ts, int *ttot); extern void clm_finalize(); @@ -89,6 +91,7 @@ GLOBAL int clmupdate_swc; GLOBAL int clmupdate_T; GLOBAL int clmupdate_texture; GLOBAL int clmprint_swc; +GLOBAL int clmprint_da_hist_file; GLOBAL int clmprint_et; GLOBAL int clmstatevec_allcol; GLOBAL int clmstatevec_colmean; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index e654b4190..4a836e85e 100755 --- a/interface/model/common/read_enkfpar.c +++ b/interface/model/common/read_enkfpar.c @@ -80,6 +80,7 @@ void read_enkfpar(char *parname) clmupdate_T = iniparser_getint(pardict,"CLM:update_T",0); clmupdate_texture = iniparser_getint(pardict,"CLM:update_texture",0); clmprint_swc = iniparser_getint(pardict,"CLM:print_swc",0); + clmprint_da_hist_file = iniparser_getint(pardict,"CLM:print_da_hist_file",0); clmprint_et = iniparser_getint(pardict,"CLM:print_et",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); clmstatevec_colmean = iniparser_getint(pardict,"CLM:statevec_colmean",0); diff --git a/interface/model/eclm/enkf_clm_5.F90 b/interface/model/eclm/enkf_clm_5.F90 index 3ae9b829f..9e2c5e593 100644 --- a/interface/model/eclm/enkf_clm_5.F90 +++ b/interface/model/eclm/enkf_clm_5.F90 @@ -71,6 +71,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") use enkf_clm_mod, only: COMM_model_clm #if defined CLMSA use enkf_clm_mod, only: define_clm_statevec + use histFileMod, only: hist_init_da_tape #endif !!<< TSMP PDAF addition end @@ -186,6 +187,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") #if defined CLMSA call define_clm_statevec(mype) + call hist_init_da_tape() #endif diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index c9f72f47a..cccf8e410 100755 --- a/interface/model/eclm/enkf_clm_mod_5.F90 +++ b/interface/model/eclm/enkf_clm_mod_5.F90 @@ -56,6 +56,7 @@ module enkf_clm_mod integer(c_int),bind(C,name="clmupdate_T") :: clmupdate_T ! by hcp integer(c_int),bind(C,name="clmupdate_texture") :: clmupdate_texture integer(c_int),bind(C,name="clmprint_swc") :: clmprint_swc + integer(c_int),bind(C,name="clmprint_da_hist_file") :: clmprint_da_hist_file #endif integer(c_int),bind(C,name="clmprint_et") :: clmprint_et integer(c_int),bind(C,name="clmstatevec_allcol") :: clmstatevec_allcol diff --git a/interface/model/eclm/hist_write_pdaf.F90 b/interface/model/eclm/hist_write_pdaf.F90 new file mode 100644 index 000000000..b743423fc --- /dev/null +++ b/interface/model/eclm/hist_write_pdaf.F90 @@ -0,0 +1,86 @@ +!------------------------------------------------------------------------------------------- +!Copyright (c) 2026 by HPSCTerrSys (Forschungszentrum Juelich GmbH) +! +!This file is part of TSMP-PDAF +! +!TSMP-PDAF is free software: you can redistribute it and/or modify +!it under the terms of the GNU Lesser General Public License as published by +!the Free Software Foundation, either version 3 of the License, or +!(at your option) any later version. +! +!TSMP-PDAF is distributed in the hope that it will be useful, +!but WITHOUT ANY WARRANTY; without even the implied warranty of +!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +!GNU Lesser General Public License for more details. +! +!You should have received a copy of the GNU Lesser General Public License +!along with TSMP-PDAF. If not, see . +!------------------------------------------------------------------------------------------- +! +! +!------------------------------------------------------------------------------------------- +! hist_write_pdaf.F90: Write CLM history file snapshots around a PDAF DA update. +! +! Two C-callable entry points: +! clm_hist_write_da_before -- call before update_clm(); captures forecast state +! clm_hist_write_da_after -- call after update_clm(); captures analysis state +! +! Both set the da_tape_phase label (used in the filename), populate the dedicated +! DA tape buffers with the current model state, and trigger hist_htapes_wrapup +! with da_call=.true. so that only the DA tape is written. +! Each call produces one time-stamped netCDF file (mfilt=1): +! caseid.clm2.da_bef.YYYY-MM-DD-SSSSS.nc +! caseid.clm2.da_aft.YYYY-MM-DD-SSSSS.nc +!------------------------------------------------------------------------------------------- + +#if defined CLMSA + +! --------------------------------------------------------------------------- +! Internal helper: set phase, update DA tape buffer, and write. +! --------------------------------------------------------------------------- +subroutine clm_hist_write_da(phase) + + use decompMod, only : get_proc_bounds, bounds_type + use clm_varpar, only : nlevgrnd + use clm_instMod, only : soilstate_inst + use histFileMod, only : hist_set_da_tape_phase, hist_update_hbuf_da, & + hist_htapes_wrapup + + implicit none + + character(len=*), intent(in) :: phase + + type(bounds_type) :: bounds + + call get_proc_bounds(bounds) + call hist_set_da_tape_phase(phase) + call hist_update_hbuf_da(bounds) + call hist_htapes_wrapup( & + rstwr = .false., & + nlend = .false., & + bounds = bounds, & + watsat_col = soilstate_inst%watsat_col(bounds%begc:bounds%endc, 1:nlevgrnd), & + sucsat_col = soilstate_inst%sucsat_col(bounds%begc:bounds%endc, 1:nlevgrnd), & + bsw_col = soilstate_inst%bsw_col(bounds%begc:bounds%endc, 1:nlevgrnd), & + hksat_col = soilstate_inst%hksat_col(bounds%begc:bounds%endc, 1:nlevgrnd), & + da_call = .true.) + +end subroutine clm_hist_write_da + +! --------------------------------------------------------------------------- +! C-callable: snapshot before the DA update (forecast / prior state) +! --------------------------------------------------------------------------- +subroutine clm_hist_write_da_before() bind(C, name="clm_hist_write_da_before") + implicit none + call clm_hist_write_da('bef') +end subroutine clm_hist_write_da_before + +! --------------------------------------------------------------------------- +! C-callable: snapshot after the DA update (analysis / posterior state) +! --------------------------------------------------------------------------- +subroutine clm_hist_write_da_after() bind(C, name="clm_hist_write_da_after") + implicit none + call clm_hist_write_da('aft') +end subroutine clm_hist_write_da_after + +#endif diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 56687e032..5f1cae9a8 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -199,7 +199,13 @@ void update_tsmp(){ #if defined CLMSA if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ + if(clmprint_da_hist_file == 1){ + clm_hist_write_da_before(); /* CLM history snapshot before DA update */ + } update_clm(&tstartcycle, &mype_world); + if(clmprint_da_hist_file == 1){ + clm_hist_write_da_after(); /* CLM history snapshot after DA update */ + } if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ print_update_clm(&tcycle, &total_steps); } From ae554d6513792cc1f95e932f8243e787143c55ef Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Mon, 27 Apr 2026 22:15:19 +0200 Subject: [PATCH 2/3] bugfix: DA history file only for eCLM fixing CLM3.5-PDAF undefined reference error --- interface/model/common/enkf.h | 2 ++ interface/model/wrapper_tsmp.c | 4 ++++ 2 files changed, 6 insertions(+) diff --git a/interface/model/common/enkf.h b/interface/model/common/enkf.h index 4b77cd840..bf8e36dda 100755 --- a/interface/model/common/enkf.h +++ b/interface/model/common/enkf.h @@ -43,9 +43,11 @@ extern void clm_advance(int *ntstep, int *tstartcycle, int *mype); extern void update_clm(int *tstartcycle, int *mype); #if defined CLMSA extern void print_update_clm(int *ts, int *ttot); +#if defined CLMFIVE extern void clm_hist_write_da_before(); extern void clm_hist_write_da_after(); #endif +#endif extern void write_clm_statistics(int *ts, int *ttot); extern void clm_finalize(); extern void cosmo_init(int pdaf_id); diff --git a/interface/model/wrapper_tsmp.c b/interface/model/wrapper_tsmp.c index 5f1cae9a8..d67022968 100644 --- a/interface/model/wrapper_tsmp.c +++ b/interface/model/wrapper_tsmp.c @@ -199,13 +199,17 @@ void update_tsmp(){ #if defined CLMSA if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0))){ +#if defined CLMFIVE if(clmprint_da_hist_file == 1){ clm_hist_write_da_before(); /* CLM history snapshot before DA update */ } +#endif update_clm(&tstartcycle, &mype_world); +#if defined CLMFIVE if(clmprint_da_hist_file == 1){ clm_hist_write_da_after(); /* CLM history snapshot after DA update */ } +#endif if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){ print_update_clm(&tcycle, &total_steps); } From 1ae719f0b6e81b186bde15e1f509df78c6d8b8fa Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 30 Apr 2026 17:24:30 +0200 Subject: [PATCH 3/3] fix: initialize the da_tape only, when CLM:print_da_hist_file is set --- interface/model/eclm/enkf_clm_5.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/interface/model/eclm/enkf_clm_5.F90 b/interface/model/eclm/enkf_clm_5.F90 index 54311f524..e2a460343 100644 --- a/interface/model/eclm/enkf_clm_5.F90 +++ b/interface/model/eclm/enkf_clm_5.F90 @@ -70,7 +70,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") use, intrinsic :: iso_C_binding, only: c_char, c_int use enkf_clm_mod, only: COMM_model_clm #if defined CLMSA - use enkf_clm_mod, only: define_clm_statevec + use enkf_clm_mod, only: define_clm_statevec, clmprint_da_hist_file use histFileMod, only: hist_init_da_tape #endif use clm_varcon, only: averaging_var @@ -189,7 +189,7 @@ subroutine clm_init(finname, pdaf_id, pdaf_max, mype) bind(C,name="clm_init") #if defined CLMSA averaging_var=0 call define_clm_statevec(mype) - call hist_init_da_tape() + if (clmprint_da_hist_file == 1) call hist_init_da_tape() #endif