diff --git a/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 b/bldsva/intf_DA/pdaf/framework/mod_read_obs.F90 index 4908d1f9..63cc7069 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/mod_tsmp.F90 b/bldsva/intf_DA/pdaf/framework/mod_tsmp.F90 index c80a22f4..311444ef 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,8 @@ 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) :: 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 a643680a..c08ba4c9 100644 --- a/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 +++ b/bldsva/intf_DA/pdaf/framework/next_observation_pdaf.F90 @@ -57,10 +57,14 @@ 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: da_interval_final + USE mod_tsmp, ONLY: flexible_da_interval USE mod_assimilation, & ONLY: obs_filename use mod_read_obs, & only: check_n_observationfile + use mod_read_obs, ONLY: check_n_observationfile_da_interval IMPLICIT NONE ! !ARGUMENTS: @@ -78,6 +82,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 @@ -100,14 +106,76 @@ SUBROUTINE next_observation_pdaf(stepnow, nsteps, doexit, time) do !nsteps = nsteps + delt_obs counter = counter + delt_obs + + ! Exit if past last observation file !if(counter>total_steps) exit - if(counter>(total_steps+toffset)) exit + if(counter>(total_steps+toffset)) then + exit + end if + + ! Check observation file #counter for observations write(fn, '(a, i5.5)') trim(obs_filename)//'.', counter call check_n_observationfile(fn,no_obs) - if(no_obs>0) exit + + ! Exit loop if observation file contains observations + if(no_obs>0) then + exit + end if + end do + + ! Set number of steps for PDAF nsteps = counter - stepnow + ! flexible_da_interval should be input (0/1) + if(flexible_da_interval.eq.1) then + +#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 + + 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) + end if + +#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 + + 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 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 ce5dd918..d3eb8a29 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; @@ -130,6 +131,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 e391f916..8e3894c7 100755 --- a/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c +++ b/bldsva/intf_DA/pdaf/model/common/read_enkfpar.c @@ -96,7 +96,11 @@ 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); 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); @@ -105,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); } } @@ -133,6 +160,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); diff --git a/doc/content/setup_tsmp/input_cmd.md b/doc/content/setup_tsmp/input_cmd.md index 3d307865..f8a6a7b4 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 987b15de..3e848dc8 100644 --- a/doc/content/setup_tsmp/input_enkfpf.md +++ b/doc/content/setup_tsmp/input_enkfpf.md @@ -70,7 +70,11 @@ dtmult = outdir = "" nreal = startreal = +total_steps = +tstartcycle = da_interval = +da_interval_final = +flexible_da_interval = stat_dumpoffset = point_obs = crns_flag = @@ -97,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 ### @@ -165,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}} @@ -617,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, @@ -635,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 @@ -670,6 +709,35 @@ 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. + +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 +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 @@ -815,6 +883,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 688968c0..c7ded102 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