From 7e21d4920499eb9b1f6a3ab44881b01fa1575415 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 11:31:47 +0200 Subject: [PATCH 01/10] next_observation_pdaf: more doc in loop --- bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index a643680af..5ca6fbdae 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -100,10 +100,16 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) do !nsteps = nsteps + delt_obs counter = counter + delt_obs + + ! Exit if at end !if(counter>total_steps) exit if(counter>(total_steps+toffset)) exit + + ! Check id observation file of counter contains observations write(fn, '(a, i5.5)') trim(obs_filename)//'.', counter call check_n_observationfile(fn,no_obs) + + ! Exit loop if observation file contains observations if(no_obs>0) exit end do nsteps = counter - stepnow From 3766ea835745a3a4d06ce426635adf1ac67e9abe Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 11:33:05 +0200 Subject: [PATCH 02/10] next_observation_pdaf: long if-statements since more code to come --- bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index 5ca6fbdae..fb47c23e9 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -103,14 +103,19 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! Exit if at end !if(counter>total_steps) exit - if(counter>(total_steps+toffset)) exit + if(counter>(total_steps+toffset)) then + exit + end if ! Check id observation file of counter contains observations write(fn, '(a, i5.5)') trim(obs_filename)//'.', counter call check_n_observationfile(fn,no_obs) ! Exit loop if observation file contains observations - if(no_obs>0) exit + if(no_obs>0) then + exit + end if + end do nsteps = counter - stepnow From c46701d6b3e78d7ae9124d8094a85101d1c3a15b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 11:50:57 +0200 Subject: [PATCH 03/10] next_observation_pdaf: adapt doc-comments in loop --- bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index fb47c23e9..591b2e685 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -101,13 +101,13 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) !nsteps = nsteps + delt_obs counter = counter + delt_obs - ! Exit if at end + ! Exit if past last observation file !if(counter>total_steps) exit if(counter>(total_steps+toffset)) then exit end if - ! Check id observation file of counter contains observations + ! Check observation file #counter for observations write(fn, '(a, i5.5)') trim(obs_filename)//'.', counter call check_n_observationfile(fn,no_obs) @@ -117,6 +117,8 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) end if end do + + ! Set number of steps for PDAF nsteps = counter - stepnow if (mype_world==0 .and. screen > 2) then From 0a1dda7cf5cc9ab326ff013321c183c831afd061 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 12:13:45 +0200 Subject: [PATCH 04/10] next_observation_pdaf: first implementation of da_interval_update so far only a dummy update --- bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 | 2 + .../pdaf/framework/next_observation_pdaf.F90 | 39 +++++++++++++++++++ bldsva/intf_DA/pdaf/model/common/enkf.h | 1 + .../intf_DA/pdaf/model/common/read_enkfpar.c | 20 ++++++++++ 4 files changed, 62 insertions(+) diff --git a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 index c80a22f4e..27bf94b06 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 +++ b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 @@ -26,6 +26,7 @@ module mod_tsmp use iso_c_binding integer(c_int) , bind(c) :: enkf_subvecsize, pf_statevecsize, nprocpf, nprocclm, nproccosmo + integer(c_int) , bind(c) :: flexible_da_interval integer(c_int) , bind(c) :: point_obs integer(c_int) , bind(c) :: is_dampfac_state_time_dependent integer(c_int) , bind(c) :: is_dampfac_param_time_dependent @@ -49,6 +50,7 @@ module mod_tsmp integer(c_int), pointer :: idx_map_subvec2state_fortran(:) type(c_ptr), bind(c) :: soilay real(c_double), pointer :: soilay_fortran(:) + real(c_double),bind(C) :: da_interval real(c_double),bind(C) :: dampfac_state_time_dependent real(c_double),bind(C) :: dampfac_param_time_dependent diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index 591b2e685..d14135e09 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -57,6 +57,8 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ONLY: mype_world USE mod_tsmp, & ONLY: total_steps + USE mod_tsmp, ONLY: da_interval + USE mod_tsmp, ONLY: flexible_da_interval USE mod_assimilation, & ONLY: obs_filename use mod_read_obs, & @@ -78,6 +80,8 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) integer :: no_obs=0 character (len = 110) :: fn !kuw end + + REAL :: da_interval_new time = 0.0 ! Not used in fully-parallel implementation variant doexit = 0 @@ -121,6 +125,41 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! Set number of steps for PDAF nsteps = counter - stepnow + ! flexible_da_interval should be input (0/1) + if(flexible_da_interval.eq.1) then + + ! total_steps can be input via `PF:simtime` + + ! Input: da_interval_final + + ! Todo: Add da_interval, da_interval_final to mod_tsmp + + ! Error Check: delt_obs must be one (check in read_enkfpar) + + ! Warning Check: no_obs, should be non-zero for every observation file + + ! Initialize da_interval_variable as zero + da_interval_new = 0 + + if(counter>(total_steps+toffset)) then + ! Add final da_interval read from input + da_interval_new = 1 !da_interval_final + else + ! Set da_interval from observation files + !call check_n_observationfile_da_interval(fn,da_interval_new) + da_interval_new = 1 + end if + + ! Warning Check: nsteps should be one + + ! Error Check: that da_interval_new is not zero anymore + + ! Update da_interval + ! TODO: Add da_interval to mod_tsmp + da_interval = da_interval_new + + end if + if (mype_world==0 .and. screen > 2) then write(*,*)'TSMP-PDAF (next_observation_pdaf.F90) stepnow: ',stepnow write(*,*)'TSMP-PDAF (next_observation_pdaf.F90) no_obs, nsteps, counter: ',no_obs,nsteps,counter diff --git a/bldsva/intf_DA/pdaf/model/common/enkf.h b/bldsva/intf_DA/pdaf/model/common/enkf.h index fba54ce56..9950d3405 100755 --- a/bldsva/intf_DA/pdaf/model/common/enkf.h +++ b/bldsva/intf_DA/pdaf/model/common/enkf.h @@ -66,6 +66,7 @@ GLOBAL int nprocclm; GLOBAL int nproccosmo; GLOBAL int nreal; GLOBAL int startreal; +GLOBAL int flexible_da_interval; GLOBAL int total_steps; GLOBAL int tcycle; GLOBAL int tstartcycle; diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c index dc8046a55..cd5133a5f 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c @@ -96,6 +96,7 @@ void read_enkfpar(char *parname) nreal = iniparser_getint(pardict,"DA:nreal",0); startreal = iniparser_getint(pardict,"DA:startreal",0); da_interval = iniparser_getdouble(pardict,"DA:da_interval",1); + flexible_da_interval = iniparser_getint(pardict,"DA:flexible_da_interval",0); stat_dumpoffset = iniparser_getint(pardict,"DA:stat_dumpoffset",0); screen_wrapper = iniparser_getint(pardict,"DA:screen_wrapper",1); point_obs = iniparser_getint(pardict,"DA:point_obs",1); @@ -132,6 +133,25 @@ void read_enkfpar(char *parname) exit(1); } + /* Check: `flexible_da_interval` must be equal to either 0 or 1 */ + /* 0: fixed da_interval (default) */ + /* 1: flexible da_interval from observation files */ + if (flexible_da_interval != 0 && flexible_da_interval != 1){ + printf("flexible_da_interval=%d\n", flexible_da_interval); + printf("Error: flexible_da_interval must be equal to either 0 or 1.\n"); + exit(1); + } + + /* Check: `da_interval` must be 1 if `flexible_da_interval` is switched on. */ + /* This way `PF:simtime` is direct input of `total_steps` + /* and `PF:starttime` is direct input of `tstartcycle`. */ + if (flexible_da_interval == 1 && da_interval != 1){ + printf("flexible_da_interval=%d\n", flexible_da_interval); + printf("da_interval=%lf\n", da_interval); + printf("Error: da_interval must be equal to 1 if flexible_da_interval is switched on.\n"); + exit(1); + } + /* Check: `npes_model = nprocpf + nprocclm + npproccosmo */ if (nprocpf + nprocclm + nproccosmo != npes_model){ printf("nprocpf=%d\n", nprocpf); From cc7d11a69ac0e639bfb238d5dfd9b6d2c67bc302 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 12:18:00 +0200 Subject: [PATCH 05/10] next_observation_pdaf: Add `da_interval_final` Final da_interval read from EnKF-input file --- bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 | 1 + .../pdaf/framework/next_observation_pdaf.F90 | 14 ++++---------- bldsva/intf_DA/pdaf/model/common/enkf.h | 1 + bldsva/intf_DA/pdaf/model/common/read_enkfpar.c | 1 + 4 files changed, 7 insertions(+), 10 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 index 27bf94b06..311444eff 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 +++ b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 @@ -51,6 +51,7 @@ module mod_tsmp type(c_ptr), bind(c) :: soilay real(c_double), pointer :: soilay_fortran(:) real(c_double),bind(C) :: da_interval + real(c_double),bind(C) :: da_interval_final real(c_double),bind(C) :: dampfac_state_time_dependent real(c_double),bind(C) :: dampfac_param_time_dependent diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index d14135e09..13009ea78 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -58,6 +58,7 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) USE mod_tsmp, & ONLY: total_steps USE mod_tsmp, ONLY: da_interval + USE mod_tsmp, ONLY: da_interval_final USE mod_tsmp, ONLY: flexible_da_interval USE mod_assimilation, & ONLY: obs_filename @@ -130,32 +131,25 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! total_steps can be input via `PF:simtime` - ! Input: da_interval_final - - ! Todo: Add da_interval, da_interval_final to mod_tsmp - ! Error Check: delt_obs must be one (check in read_enkfpar) - ! Warning Check: no_obs, should be non-zero for every observation file + ! Warning Check: nsteps should be one ! Initialize da_interval_variable as zero da_interval_new = 0 if(counter>(total_steps+toffset)) then - ! Add final da_interval read from input - da_interval_new = 1 !da_interval_final + ! Set da_interval to da_interval_final from EnKF input file + da_interval_new = da_interval_final else ! Set da_interval from observation files !call check_n_observationfile_da_interval(fn,da_interval_new) da_interval_new = 1 end if - ! Warning Check: nsteps should be one - ! Error Check: that da_interval_new is not zero anymore ! Update da_interval - ! TODO: Add da_interval to mod_tsmp da_interval = da_interval_new end if diff --git a/bldsva/intf_DA/pdaf/model/common/enkf.h b/bldsva/intf_DA/pdaf/model/common/enkf.h index 9950d3405..beff4bdc1 100755 --- a/bldsva/intf_DA/pdaf/model/common/enkf.h +++ b/bldsva/intf_DA/pdaf/model/common/enkf.h @@ -130,6 +130,7 @@ GLOBAL double dt; GLOBAL double pf_aniso_perm_y; GLOBAL double pf_aniso_perm_z; GLOBAL double da_interval; +GLOBAL double da_interval_final; GLOBAL double pf_dampfac_param; GLOBAL double pf_dampfac_state; GLOBAL double dampfac_state_time_dependent; diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c index cd5133a5f..baad54895 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c @@ -96,6 +96,7 @@ void read_enkfpar(char *parname) nreal = iniparser_getint(pardict,"DA:nreal",0); startreal = iniparser_getint(pardict,"DA:startreal",0); da_interval = iniparser_getdouble(pardict,"DA:da_interval",1); + da_interval_final = iniparser_getdouble(pardict,"DA:da_interval_final",1); flexible_da_interval = iniparser_getint(pardict,"DA:flexible_da_interval",0); stat_dumpoffset = iniparser_getint(pardict,"DA:stat_dumpoffset",0); screen_wrapper = iniparser_getint(pardict,"DA:screen_wrapper",1); From 10ab3c1a5f6797ac6d286d00b9bddd29de19191d Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 12:36:00 +0200 Subject: [PATCH 06/10] next_observation_pdaf: Add `da_interval` from observation file --- .../intf_DA/pdaf/framework/mod_read_obs.F90 | 41 +++++++++++++++++++ .../pdaf/framework/next_observation_pdaf.F90 | 15 +++++-- 2 files changed, 52 insertions(+), 4 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 b/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 index 4908d1f92..63cc70691 100755 --- a/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 +++ b/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 @@ -585,6 +585,47 @@ subroutine check_n_observationfile(fn,nn) end subroutine check_n_observationfile + !> @author Yorck Ewerdwalbesloh, Johannes Keller + !> @date 11.09.2023 + !> @brief Return data assimilation interval from file + !> @param[in] fn Filename of the observation file + !> @param[out] aa new da_interval (number of time steps until next assimilation time step) + !> @details + !> Reads the content of the variable name `da_interval` from NetCDF + !> file `fn` using subroutines from the NetCDF module. + !> The result is returned in `aa`. + !> + !> The result is used to adapt the da_interval until the next observation file. + !> + !> Adapted for TSMP-PDAF by Johannes Keller, 28.05.2025. + subroutine check_n_observationfile_da_interval(fn,aa) + use shr_kind_mod, only: r8 => shr_kind_r8 + use netcdf, only: nf90_max_name, nf90_open, nf90_nowrite, & + nf90_inq_varid, nf90_get_var, nf90_close, nf90_noerr + + implicit none + + character(len=*),intent(in) :: fn + real, intent(out) :: aa + + + integer :: ncid, varid, status !,dimid + character (len = *), parameter :: varname = "da_interval" + real(r8) :: dtime ! land model time step (sec) + + !character (len = *), parameter :: dim_name = "dim_obs" + !character(len = nf90_max_name) :: recorddimname + + call check(nf90_open(fn, nf90_nowrite, ncid)) + !call check(nf90_inq_dimid(ncid, dim_name, dimid)) + !call check(nf90_inquire_dimension(ncid, dimid, recorddimname, nn)) + + call check( nf90_inq_varid(ncid, varname, varid)) + call check( nf90_get_var(ncid, varid, aa) ) + call check(nf90_close(ncid)) + + end subroutine check_n_observationfile_da_interval + !> @author Wolfgang Kurtz, Guowei He, Mukund Pondkule !> @date 03.03.2023 !> @brief Error handling for netCDF commands diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index 13009ea78..87e2f069f 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -64,6 +64,7 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ONLY: obs_filename use mod_read_obs, & only: check_n_observationfile + use mod_read_obs, ONLY: check_n_observationfile_da_interval IMPLICIT NONE ! !ARGUMENTS: @@ -136,18 +137,24 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! Warning Check: nsteps should be one ! Initialize da_interval_variable as zero - da_interval_new = 0 + da_interval_new = 0.0 if(counter>(total_steps+toffset)) then ! Set da_interval to da_interval_final from EnKF input file da_interval_new = da_interval_final else ! Set da_interval from observation files - !call check_n_observationfile_da_interval(fn,da_interval_new) - da_interval_new = 1 + call check_n_observationfile_da_interval(fn,da_interval_new) end if - ! Error Check: that da_interval_new is not zero anymore +#ifdef PDAF_DEBUG + ! Error Check: da_interval_new should be set to at least one + if(da_interval_new < 1.0) then + write(*,'(a,es22.15)') "da_interval_new = ", da_interval_new + write(*,'(a)') "da_interval_new is too small, should be minimum of one" + stop "Stopped from incorrect da_interval_new" + end if +#endif ! Update da_interval da_interval = da_interval_new From 3a166f93ff39a27431072cdbda43e3f8945ca9b3 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 12:46:13 +0200 Subject: [PATCH 07/10] next_observation_pdaf: Checks for flexible time stepping --- .../pdaf/framework/next_observation_pdaf.F90 | 31 +++++++++++++++++-- 1 file changed, 28 insertions(+), 3 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index 87e2f069f..261ad9655 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -130,11 +130,32 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! flexible_da_interval should be input (0/1) if(flexible_da_interval.eq.1) then - ! total_steps can be input via `PF:simtime` + ! Document: EnKF-Parameter-File Input `flexible_da_interval` - ! Error Check: delt_obs must be one (check in read_enkfpar) + ! Document: total_steps can be input via `PF:simtime` - ! Warning Check: nsteps should be one + ! Document: Input da_interval_final + + ! Document: Input restriction for delt_obs + + ! Document: Observation file input `da_interval` + +#ifdef PDAF_DEBUG + ! Error Check: delt_obs must be one for flexible time stepping + if (delt_obs /= 1) then + write(*,'(a,i10)') "delt_obs = ", delt_obs + write(*,'(a)') "delt_obs must be one for flexible time stepping" + stop "Stopped from incorrect delt_obs" + end if + + ! Warning: nsteps should be one + if(nsteps > 1) then + write(*,'(a,i10)') "WARNING: nsteps = ", nsteps + write(*,'(a)') "WARNING: nsteps should be one for flexible time stepping" + write(*,'(a)') "WARNING: Any time differences can be encoded in observation files" + write(*,'(a)') "WARNING: using the variable da_interval." + end if +#endif ! Initialize da_interval_variable as zero da_interval_new = 0.0 @@ -159,6 +180,10 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! Update da_interval da_interval = da_interval_new + if (mype_world==0 .and. screen > 2) then + write(*,'(a,es22.15)')'TSMP-PDAF (next_observation_pdaf.F90) da_interval: ', da_interval + end if + end if if (mype_world==0 .and. screen > 2) then From d8f336cf946e54956f0c813eee3dfd6df2ef065a Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Wed, 28 May 2025 13:21:33 +0200 Subject: [PATCH 08/10] documentation: flexible time stepping --- .../pdaf/framework/next_observation_pdaf.F90 | 10 ------- doc/content/setup_tsmp/input_cmd.md | 2 ++ doc/content/setup_tsmp/input_enkfpf.md | 28 +++++++++++++++++++ doc/content/setup_tsmp/input_obs.md | 5 ++++ 4 files changed, 35 insertions(+), 10 deletions(-) diff --git a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 index 261ad9655..c08ba4c9c 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -130,16 +130,6 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) ! flexible_da_interval should be input (0/1) if(flexible_da_interval.eq.1) then - ! Document: EnKF-Parameter-File Input `flexible_da_interval` - - ! Document: total_steps can be input via `PF:simtime` - - ! Document: Input da_interval_final - - ! Document: Input restriction for delt_obs - - ! Document: Observation file input `da_interval` - #ifdef PDAF_DEBUG ! Error Check: delt_obs must be one for flexible time stepping if (delt_obs /= 1) then diff --git a/doc/content/setup_tsmp/input_cmd.md b/doc/content/setup_tsmp/input_cmd.md index 3d3078654..f8a6a7b4b 100644 --- a/doc/content/setup_tsmp/input_cmd.md +++ b/doc/content/setup_tsmp/input_cmd.md @@ -73,6 +73,8 @@ In general, `delt_obs` should be as small as possible, in order to avoid performance loss. See remark in [`[DA]da_interval`](./input_enkfpf.md#dada_interval). +For flexible time stepping, `delt_obs` must be one. + ## screen ## `screen` (integer) Control verbosity of PDAF diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index 4d1082a1a..2ce59bfe4 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -70,6 +70,8 @@ outdir = "" nreal = startreal = da_interval = +da_interval_final = +flexible_da_interval = stat_dumpoffset = point_obs = crns_flag = @@ -653,6 +655,30 @@ called, assembling EnKF state vectors and calling the PDAF library. Maximizing `da_interval`, minimizes the number of these calls and thus reduces compute time. +### DA:da_interval_final ### + +`DA:da_interval_final`: (float) Final `da_interval` for flexible time +stepping using observation files. + +Used only if `DA:flexible_da_interval` is switched on. + +### DA:flexible_da_interval ### + +`DA:flexible_da_interval`: (integer) Switch for flexible time stepping +using observation files. + +- `0`: default, fixed time stepping, single `da_interval` +- `1`: flexible time stepping, `da_interval` read from next + observation file + +- Set `total_steps` using `PF:simtime` (`DA:da_interval` must be `1`) +- Command line input [`delt_obs`](./input_cmd.md#delt_obs) must be + `1`! + +In observation files, the new input variable +[`da_interval`](./input_obs.md#da_interval) determines the +`da_interval` before the observation file in question. + ### DA:stat_dumpoffset ### `DA:stat_dumpoffset`: File number offset for the data assimilation @@ -797,6 +823,8 @@ Default: 0, output turned off. | | `nreal` | 0 | | | `outdir` | \- | | | `da_interval` | 1 | + | | `da_interval_final` | 1 | + | | `flexible_da_interval` | 0 | | | `stat_dumpoffset` | 0 | | | `screen_wrapper` | 1 | | | `point_obs` | 1 | diff --git a/doc/content/setup_tsmp/input_obs.md b/doc/content/setup_tsmp/input_obs.md index 688968c00..c7ded1028 100644 --- a/doc/content/setup_tsmp/input_obs.md +++ b/doc/content/setup_tsmp/input_obs.md @@ -393,6 +393,11 @@ f.variables["no_obs"].assignValue(0) f.close() ``` +### da_interval ### + +`da_interval`: (float) Value for +[`da_interval`](./input_enkfpf.md#dada_interval) used in the +assimilation cycle leading up to the observation file. ## Specifying type of observation at compile time ## The following preprocessor variables let the user specify the expected From 595912a9e5614124707ea0c170afb39afdc07a8b Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Jun 2025 15:33:14 +0200 Subject: [PATCH 09/10] doc: note about `da_interval_final` --- doc/content/setup_tsmp/input_enkfpf.md | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index a05728a68..0250d44b9 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -679,6 +679,11 @@ stepping using observation files. Used only if `DA:flexible_da_interval` is switched on. +Note: If the last observation file contains observations, +`da_interval_final` is set, but never used, as the simulation +currently finishes exactly at the time the observations from the last +observation file is assimilated. + ### DA:flexible_da_interval ### `DA:flexible_da_interval`: (integer) Switch for flexible time stepping From 6386a988d1f42505056fe3f1cc5bb50cba7b3d97 Mon Sep 17 00:00:00 2001 From: Johannes Keller Date: Thu, 5 Jun 2025 16:50:54 +0200 Subject: [PATCH 10/10] New inputs `DA:total_steps` and `DA:tstartcycle` Option to set `total_steps` and `tstartcycle` directly. Before / by default, they are set through `PF:simtime` and `PF:starttime`, respectively. --- .../intf_DA/pdaf/model/common/read_enkfpar.c | 33 ++++++++++-- doc/content/setup_tsmp/input_enkfpf.md | 53 ++++++++++++++++--- 2 files changed, 74 insertions(+), 12 deletions(-) diff --git a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c index f3c996d0a..8e3894c79 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c @@ -96,6 +96,8 @@ void read_enkfpar(char *parname) strcpy(outdir,string); nreal = iniparser_getint(pardict,"DA:nreal",0); startreal = iniparser_getint(pardict,"DA:startreal",0); + total_steps = iniparser_getint(pardict,"DA:total_steps",-1); + tstartcycle = iniparser_getint(pardict,"DA:tstartcycle",-1); da_interval = iniparser_getdouble(pardict,"DA:da_interval",1); da_interval_final = iniparser_getdouble(pardict,"DA:da_interval_final",1); flexible_da_interval = iniparser_getint(pardict,"DA:flexible_da_interval",0); @@ -107,15 +109,38 @@ void read_enkfpar(char *parname) da_crns_depth_tol = iniparser_getdouble(pardict,"DA:da_crns_depth_tol",0.01); clmcrns_bd = iniparser_getdouble(pardict, "DA:crns_bd", -1.0); da_print_obs_index = iniparser_getint(pardict,"DA:print_obs_index",0); - total_steps = (int) (t_sim/da_interval); - tstartcycle = (int) (t_start/da_interval); - /* print inputs / debug output for data assimilation settings */ + /* Debug output */ if (mype_world == 0) { if (screen_wrapper > 0) { printf("TSMP-PDAF-WRAPPER mype(w)=%5d read_enkfpar: [DA]\n",mype_world); printf("TSMP-PDAF-WRAPPER mype(w)=%5d ------------------\n",mype_world); - printf("TSMP-PDAF-WRAPPER mype(w)=%5d t_sim = %lf | da_interval = %lf | total_steps = %d\n",mype_world,t_sim,da_interval,total_steps); + } + } + + /* Set total_steps from PF:simtime if it has not been set + directly */ + if (total_steps == -1) { + total_steps = (int) (t_sim/da_interval); + + /* Debug output */ + if (mype_world == 0) { + if (screen_wrapper > 0) { + printf("TSMP-PDAF-WRAPPER mype(w)=%5d total_steps set based on t_sim = %lf\n",mype_world,t_sim); + } + } + } + + /* Set tstartcycle from PF:starttime if it has not been set + directly */ + if (tstartcycle == -1) { + tstartcycle = (int) (t_start/da_interval); + } + + /* Debug output */ + if (mype_world == 0) { + if (screen_wrapper > 0) { + printf("TSMP-PDAF-WRAPPER mype(w)=%5d da_interval = %lf | total_steps = %d\n",mype_world,da_interval,total_steps); printf("TSMP-PDAF-WRAPPER mype(w)=%5d nreal = %d | n_modeltasks = %d\n",mype_world,nreal,n_modeltasks); } } diff --git a/doc/content/setup_tsmp/input_enkfpf.md b/doc/content/setup_tsmp/input_enkfpf.md index 0250d44b9..3e848dc85 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -70,6 +70,8 @@ dtmult = outdir = "" nreal = startreal = +total_steps = +tstartcycle = da_interval = da_interval_final = flexible_da_interval = @@ -99,8 +101,12 @@ available in a single `COMM_model`. Each `COMM_model` contains ### PF:starttime ### -`PF:starttime`: (real) ParFlow start time. Must match with the -specifications in the `*.pfidb` input (`TimingInfo.StartTime`). +`PF:starttime`: (real) ParFlow start time. + +Not used if `DA:tstartcycle` is set. + +If used: Must match with the specifications in the `*.pfidb` input +(`TimingInfo.StartTime`). ### PF:dt ### @@ -167,11 +173,14 @@ Deprecated. Sets `PF:simtime`. `PF:simtime`: (real) Total simulation time (in terms of ParFlow timing). -For all simulations (including CLMSA): `PF:simtime` is used in -conjunction with `DA:da_interval` to determine `total_steps`, the -total number of iterations of the main data assimilation loop. Each -iteration of the main data assimilation loop consists of one forward -simulation and one data assimilation step. +Not used if `DA:total_steps` is set. + +If `DA:total_steps` is not set, the following holds for all +simulations (including CLMSA): `PF:simtime` is used in conjunction +with `DA:da_interval` to determine `total_steps`, the total number of +iterations of the main data assimilation loop. Each iteration of the +main data assimilation loop consists of one forward simulation and one +data assimilation step. \begin{align*} \mathtt{total\_steps} &= \frac{\mathtt{PF:simtime}}{\mathtt{DA:da\_interval}} @@ -619,6 +628,32 @@ errors. Recommendation: Use at least an ensemble size of 4. `DA:startreal`: (integer) Added to suffix-numbers for input file creation. +### DA:total_steps ### + +`DA:total_steps`: (integer) Number of observation Intervals. + +Together with `DA:da_interval` and `DA:startcycle`, `DA:total_steps` +determines the simulation time. + +If not set directly, `DA:total_steps` is determined from `PF:simtime`. + +`DA:total_steps` must be in sync with timing information from +component models. + +- ParFlow: `DA:total_steps` times `DA:da_interval` should be + `TimingInfo.StopTime` minus `TimingInfo.StartTime` (from `*.pfidb`). + +- CLM: `DA:total_steps` times `DA:da_interval` should be the number of + CLM time steps in `stop_ymd` minus `start_ymd start_tod`. To be + safe, CLM's stoptime can be chosen later, such that it will be + stopped by the TSMP-PDAF stop alarm. + +### DA:tstartcycle ### + +`DA:tstartcycle`: (integer) First observation cycle (file) to use. + +This should be zero. Other values are currently not supported. + ### DA:da_interval ### `DA:da_interval`: (double) Time interval (units of ParFlow timing, @@ -637,7 +672,9 @@ assimilation. For CLM simulations, `DA:da_interval` determines the number of CLM-time-steps between two assimilations together with `PF:dt`. See -the section of `PF:dt` for more detail. +the section of `PF:dt` for more detail. In particular, for CLM +simulations and `PF:dt==1.0`, `DA:da_interval` is in units of CLM time +steps. One exception is that the observation file is empty at a data assimilation step. Then, no data assimilation takes place and the next