diff --git a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md b/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md index a3902a84..41ad67ea 100644 --- a/docs/users_guide/running_tsmp_pdaf/input_enkfpf.md +++ b/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.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: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 | diff --git a/interface/model/Makefile b/interface/model/Makefile index 185b906a..35859efa 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 7462a6f1..61506554 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 a734c335..06a1067f 100755 --- a/interface/model/common/enkf.h +++ b/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; diff --git a/interface/model/common/read_enkfpar.c b/interface/model/common/read_enkfpar.c index 27b26c46..cc4bd284 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); clmprint_inc = iniparser_getint(pardict,"CLM:print_inc",0); clmstatevec_allcol = iniparser_getint(pardict,"CLM:statevec_allcol",0); diff --git a/interface/model/eclm/enkf_clm_5.F90 b/interface/model/eclm/enkf_clm_5.F90 index dd7b43a8..e2a46034 100644 --- a/interface/model/eclm/enkf_clm_5.F90 +++ b/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 diff --git a/interface/model/eclm/enkf_clm_mod_5.F90 b/interface/model/eclm/enkf_clm_mod_5.F90 index bc2e4c59..adde3aae 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 ! Yorck integer(c_int),bind(C,name="clmupdate_tws") :: clmupdate_tws diff --git a/interface/model/eclm/hist_write_pdaf.F90 b/interface/model/eclm/hist_write_pdaf.F90 new file mode 100644 index 00000000..b743423f --- /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 87c2d093..18d615fb 100644 --- a/interface/model/wrapper_tsmp.c +++ b/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);