Flexible timing in TSMP-PDAF

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

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
 

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

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;

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);

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