Add CLM history file output around DA update

docs/users_guide/running_tsmp_pdaf/input_enkfpf.md

@@ -56,6 +56,7 @@ update_texture  =
 update_T  =
 print_swc   =
 print_et   =
+print_da_hist_file =
 print_inc   =
 statevec_allcol =
 statevec_colmean =
@@ -513,6 +514,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<inst_suffix>.da_bef.YYYY-MM-DD-SSSSS.nc` and
+`caseid.clm2<inst_suffix>.da_aft.YYYY-MM-DD-SSSSS.nc`, where the
+timestamp is the model time of the DA update. Default: `0`.
+
 ### CLM:print_inc ###
 
 `CLM:print_inc`: (integer) If set to `1`, the analysis increment
@@ -1026,6 +1038,7 @@ Default: 0, output turned off.
  |           | `update_swc`            | 1             |
  |           | `print_swc`             | 0             |
  |           | `print_et`              | 0             |
+ |           | `print_da_hist_file`    | 0             |
  |           | `print_inc`             | 0             |
  |           | `statevec_allcol`       | 0             |
  |           | `statevec_colmean`      | 0             |

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

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

interface/model/common/enkf.h

@@ -43,6 +43,10 @@ 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
 extern void print_inc_clm();
 #endif
 extern void write_clm_statistics(int *ts, int *ttot);
@@ -92,6 +96,7 @@ GLOBAL int clmupdate_T;
 GLOBAL int clmupdate_texture;
 GLOBAL int clmupdate_tws;
 GLOBAL int clmprint_swc;
+GLOBAL int clmprint_da_hist_file;
 GLOBAL int clmprint_et;
 GLOBAL int clmprint_inc;
 GLOBAL int clmstatevec_allcol;

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);
   clmprint_inc           = iniparser_getint(pardict,"CLM:print_inc",0);
   clmstatevec_allcol    = iniparser_getint(pardict,"CLM:statevec_allcol",0);

interface/model/eclm/enkf_clm_5.F90

@@ -70,7 +70,8 @@ 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
 !!<< TSMP PDAF addition end
@@ -188,6 +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)
+  if (clmprint_da_hist_file == 1) call hist_init_da_tape()
 #endif
 
 

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
 
   ! Yorck
   integer(c_int),bind(C,name="clmupdate_tws") :: clmupdate_tws

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 <http://www.gnu.org/licenses/>.
+!-------------------------------------------------------------------------------------------
+!
+!
+!-------------------------------------------------------------------------------------------
+! 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<inst_suffix>.da_bef.YYYY-MM-DD-SSSSS.nc
+!   caseid.clm2<inst_suffix>.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

interface/model/wrapper_tsmp.c

@@ -199,8 +199,19 @@ void update_tsmp(){
 
 #if defined CLMSA
   if((model == tag_model_clm) && ((clmupdate_swc != 0) || (clmupdate_T != 0) || (clmupdate_tws != 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
 #ifndef CLMFIVE
     if(clmprint_swc == 1 || clmupdate_texture == 1 || clmupdate_texture == 2){
       print_update_clm(&tcycle, &total_steps);