Incremental Analysis Updates (IAU) with PDAF

IAU (Incremental Analysis Updates) was re-coded for PDAF V3.0 and now offers full functionality

Incremental Analysis Updates (IAU) are a method to avoid shocks to the model dynamics that could result from the increment computed in an analysis step. Instead of an instantaneous change of the model fields, IAU applies the analysis increment stepwise distributes over a number of time steps.

There are different approaches to IAU. One is what we call forward IAU. In this case the ensmeble is integrated by the model forard in time starting from the analysis time. Thus, the analysis increment is applied in the future of the analysis step. The other approach is retrospective IAU. This approach usually assimilates observation in the middle of a time window. Then, the model is restarted at the beginning of the time window and the IAU is applied over the full time window. In this approach the distributed increments are usually weighted to that they incremental changes increase toward the analysis time in the middle of the time window and decrease afterwards. An example would be for daily assimilation to compute the analysis increment in the middle of the day and then restart the ensemble integration from midnight and integrating over 24 hours.

PDAF supports both forward IAU and retrospective IAU, however, for retrospective IAU with online coupled DA the user would need to ensure that the model steps back in time consistently.

If one uses the PDAF_assimilate interface routines, one can activate IAU with a single call to PDAF_iau_init, which will be explained below. Apart from the standard settings, there are further routines that can be called for particular cases.

Activation and configuration of IAU

PDAF_iau_init

This routine initializes parameters for IAU and activates it. The routine further allocates the array in which the ensemble increments are stored. This array exists on all processes that are part of model tasks.

The routine is usually called in init_pdaf after the initialization of PDAF in PDAF_init. it is important that the routine is called by all model processes because all these processes need the information on the IAU configuration and need to allocate the increment array.

The interface is:

SUBROUTINE PDAF_iau_init(type_iau, nsteps_iau, flag)

    INTEGER, INTENT(in) :: type_iau     ! Type of IAU
                                        !   (0) no IAU
                                        !   (1) constant increment weight 1/nsteps_iau
                                        !   (2) Linear IAU weight with maximum in middle of IAU period
                                        !   (3) Zero weights for null mode (can be used to apply IAU on user side)
    INTEGER, INTENT(in) :: nsteps_iau   ! Number of time steps in IAU
    INTEGER, INTENT(out) :: flag        ! Status flag

Hints:

• type_iau=2 is linearly increasing and decreasing similar to what is used in the ocean model NEMO. This type is usually used if the IAU is applied in re-running over the previous observation period.
• type_iau=3 stores the increment information, but does not apply the increment. This permits to apply the increment in user-provided code. For this, the routine PDAF_iau_set_pointer can be used to assess the increment array within PDAF.

An example implementation can be found in tests/online_2D_parallelmodel/.

Set the increment array

PDAF_iau_init_inc

PDAF_iau_init_inc can be used by a user to fill the process-local increment array. A common use is when the IAU should be applied from the initial time of an assimilation run, for example if the increment was computed in a previous assimilation run and then stored for restarting. Since PDAF can only compute an increment in an analysis step, the user needs to provide the increment.

The routine is usually called in init_pdaf after the initialization of IAU with PDAF_iau_init. It has to be called by all processes that are model processes and one needs to provide the task-local ensemble (i.e. with local ensemble size dim_ens_l=1 for the fully parallel mode, and usually dim_ens_l>1 for the flexible parallelization mode). The routine can usually not be called in init_ens_pdaf since this routine is only executed by filter processes and not all model processes.

The interface is:

  SUBROUTINE PDAF_iau_init_inc(dim_p, dim_ens_l, ens_inc, flag)

    INTEGER, INTENT(in) :: dim_p                    ! PE-local dimension of model state
    INTEGER, INTENT(in) :: dim_ens_l                ! Task-local size of ensemble
    REAL, INTENT(in) :: ens_inc(dim_p, dim_ens_l)   ! PE-local increment ensemble
    INTEGER, INTENT(out) :: flag                    ! Status flag

Change the IAU configuration

PDAF_iau_reset

PDAF_iau_reset is used to modify the IAU type and the number of IAU time steps during a run. While `PDAF_iau_init` sets the IAU type and number of IAU steps initially, one can change these settings during an assimilation run.

The routine has to be called by all processes that are model processes. A common place is to call the routine in assimilate_pdaf after an analysis step or in distribute_state_pdaf.

The interface is:

  SUBROUTINE PDAF_iau_reset(type_iau, nsteps_iau, flag)

    INTEGER, INTENT(in) :: type_iau     ! Type of IAU
                                        !   (0) no IAU
                                        !   (1) constant increment weight 1/nsteps_iau
                                        !   (2) Linear IAU weight with maximum in middle of IAU period
                                        !   (3) Zero weights for null mode (can be used to apply IAU on user side)
    INTEGER, INTENT(in) :: nsteps_iau   ! number of time steps in IAU
    INTEGER, INTENT(out) :: flag        ! Status flag

Hints:

• It is recommended to only call this routine after or just before a forecast phase.
• Calling this routine during a forecast phase will change the IAU weights and might stop the application of IAU if nsteps_iau is below the number of time steps performed in the current forecast phase.

Set user-defined IAU weights

PDAF_iau_set_weights

This routine allows the user to provide a user-specified vector of increment weights. While PDAF_iau_init allows to choose among pre-defined weight functions, one might like to use a different function and the corresponding weights can be set here.

The routine has to be called by all processes that are model processes. A common place is to call the routine in assimilate_pdaf after an analysis step or in distribute_state_pdaf.

The interface is:

  SUBROUTINE PDAF_iau_set_weights(iweights, weights)

    INTEGER, INTENT(in) :: iweights        ! Length of weights input vector
    REAL, INTENT(in) :: weights(iweights)  ! Input weight vector

Hints:

• If iweights is different from the number of IAU steps set in PDAF_iau_init or PDAF_iau_reset only the minimum of iweights and the set IAU steps is filled with the provided weights vector.
• We recommend to only change the weights vector before a forecast phase.
• While setting different weight functions for different processes or model tasks is possible this likely leads to inconsistencies.

Set a pointer to the IAU array

PDAF_iau_set_pointer

With this routine the user can set a pointer to the PDAF-internal array of ensemble increments. This gives direct access to the increment array e.g. to analyze it or to write it into a file for restarting. The routine is different from PDAF_iau_set_inc in that here the user obtains access to the internal IAU array, while in PDAF_iau_set_inc one can only overwrite the internal array.

The routine can be called by each single process, but it only provides a pointer to the process-local part of the increment array. For domain-decomposed models, this array only includes the state vector part for the process domain. In addition, it usually only contains a sub-ensemble unless one uses the flexible parallelization mode with a single model task. For the fully parallel mode, the process(es) of a single model task only hold a single ensemble state.

The interface is:

  SUBROUTINE PDAF_iau_set_pointer(iau_ptr, flag)

    REAL, POINTER, INTENT(out) :: iau_ptr(:,:)  ! Pointer to IAU ensemble array
    INTEGER, INTENT(out)       :: flag          ! Status flag

Hints:

• In Fortran user code one has to declare
  REAL, POINTER :: iau_ptr(:,:)
  and provide this as the first argument. One does not need to allocate this pointer.

Apply the IAU increment

PDAF_iau_add_inc

This routine is only used with the flexible parallelization mode if a routine PDAF3_put_state_X (or any other PDAF*_put_state_X) routine is used. Since PDAF*_put_state_X is only called after a model state is integrated over the full forecast period, this routine cannot apply the increment. Hence, the routine PDAF_iau_add_inc has to be called at each model time step in the time stepping loop of the model to apply the increment. This routine adds the increment for the currently integrated ensemble member according to the specified IAU type and number of IAU steps.

Note: Starting from PDAF V3.0 we recommend to use the PDAF*_assimilate_X routines also for the flexible parallelization mode. Since PDAF*_assimilate_X is called at each time step, the additional call of PDAF_iau_add_inc is not required.

The routine has to be inserted in time stepping loop of the model code. It has to be called at each time step.

The interface is the following:

  SUBROUTINE PDAF_iau_add_inc(U_collect_state, U_distribute_state)

    EXTERNAL :: U_collect_state, &      ! Routine to collect a state vector
         U_distribute_state             ! Routine to distribute a state vector