wiki:AuxiliaryRoutines

Auxiliary Routines

PDAF provide a number of auxiliary routines. They provide access to PDAF-internal data, which is not available through the regular interface of PDAF.

Interface specifications

PDAF_interfaces_module.F90

This file contains the module PDAF_interfaces_module. It provides interface definitions for the routines of PDAF. If you like to use explicit Fortran interfaces, or if you have to use then in case of the smoother, you have to include the following line in the header you the calling routine:

  use PDAF_interfaces_module

Access to PDAF-internal arrays

PDAF_get_smootherens

A smoother is available for several filters (ESTKF/LESTKF/ETKF/LETKF/EnKF). This routine is called to set a Fortran pointer to the array in PDAF that stores the ensembles for smoothing. Also, the routine sets the available lag of the smoothing.

CALL PDAF_get_smootherens(sens_pointer, maxlag, status)

The variables in the interface are the following:

  • sens_pointer: The pointer to the smoother ensemble. The dimension is sens_pointer(:,:,:). Thus in the program calling PDAF_get_smootherens one has to declare REAL, POINTER :: sens_pointer(:,:,:). On output it points to the smoother ensemble.
  • maxlag: Number of lags stored in the smoother ensemble. While in the call to PDAF_init the maximum lag for the smoother is set, not all possible lags are using at the beginning of the assimilations. maxlag shows how many times were already smoothed.
  • status: Status flag. 0 for successful exit.

Important:

  • Because sens_pointer is a pointer, the call to PDAF_get_smootherens needs an explicit Fortran interface. This is provided by the Fortran module PDAF_interfaces_module. In the header part of the routine that calls PDAF_get_smootherens one has to include use PDAF_interfaces_module!
  • Using a pointer combined with an intent, i.e. using a pointer as argument, is a feature of Fortran 2003. Thus, if a too old compiler is used, it will provide an error when the routine is compiled.

Notes:

  • PDAF_get_smootherens is typically called in the prepoststep routine after the analysis step. At this time not only the filtered analysis step can be analized, but also all smoothed ensembles.
  • The first two indices of sens_pointer are identical to those in the ensemble array ens_p. Thus, the array contains state vectors in its columns. The second index is the ensemble index. The third index is the lag. Thus, if the value of the third index is fixed one can use the array sens_pointer analogous to the ensemble array ens_p.
  • For an example using PDAF_get_smootherens see the implementation for the Lorenz96 model in models/lorenz96/. The routine is called in compute_rms_smoother.F90.

PDAF_set_smootherens

A smoother is available for several filters (ESTKF/LESTKF/ETKF/LETKF/EnKF/NETF/LNETF). This routine is called to set a Fortran pointer to the array in PDAF that stores the ensembles for smoothing. In addition, it sets the available lag of the smoothing. This routine is called in the offline mode of PDAF. While in the online mode, the smoother ensemble array is filled automatically by PDAF, one has to fill it manually in the offline mode. PDAF_set_smootherens gives access to the smoother array to fill it.

CALL PDAF_set_smootherens(sens_pointer, maxlag, status)

The arguments are:

  • sens_pointer: The pointer to the smoother ensemble. The dimension is sens_pointer(:,:,:). Thus in the program calling PDAF_get_smmotherens one has to declare REAL, POINTER :: sens_pointer(:,:,:). On output it points to the smoother ensemble.
  • maxlag: Set the number of lags stored in the smoother ensemble. While in the call to PDAF_init the maximum lag for the smoother is set, not all possible lags are using at the beginning of the assimilations. maxlag says how many times were already smoothed. Both values are usually identical for the offline mode.
  • status: Status flag. 0 for successful exit.

Important:

  • Because sens_pointer is a pointer, the call to PDAF_set_smootherens needs an explicit Fortran interface. This is provided by the Fortran module PDAF_interfaces_module. In the header part of the routine that calls PDAF_set_smootherens one has to include use PDAF_interfaces_module!
  • Using a pointer combined with an intent, i.e. using a pointer as argument, is a feature of Fortran 2003. Thus, if a too old compiler is used, it will provide an error when the routine is compiled.

Notes:

  • PDAF_set_smootherens is typically called in the initialization phase of PDAF. It was to be called after the basic initialization of PDAF in PDAF_init. A possible location is to call PDAF_set_smootherens is the ensemble initialization routine U_init_ens.
  • The first two indices of sens_pointer are identical to those in the ensemble array ens_p. Thus, the array contains state vectors in its columns. The second index is the ensemble index. The third index is the lag. Thus, if the value of the third index is fixed one can use the array sens_pointer analogous to the ensemble array ens_p.

PDAF_set_ens_pointer

This routine allows a program to get a Fortran pointer to the internal ensemble array of PDAF.

CALL PDAF_set_ens_pointer(ens_pointer, status)

The arguments are:

  • ens_pointer: The pointer to the smoother ensemble. The dimension is ens_pointer(:,:). Thus in the program calling PDAF_set_ens_pointer one has to declare REAL, POINTER :: sens_pointer(:,:). On output it points to the ensemble array.
  • status: Status flag. 0 for successful exit.

Important:

  • Because ens_pointer is a pointer, the call to PDAF_set_ens_pointer needs an explicit Fortran interface. This is provided by the Fortran module PDAF_interfaces_module. In the header part of the routine that calls PDAF_set_ens_pointer one has to include use PDAF_interfaces_module!
  • Using a pointer combined with an intent, i.e. using a pointer as argument, is a feature of Fortran 2003. Thus, if a too old compiler is used, it will provide an error when the routine is compiled.

Notes:

  • PDAF_set_ens_pointer is a special routine that is never needed when the standard online or offline modes of the implementation are used. However, the routine allows to build a model that uses each column of the ensemble array to store the model fields. Thus, one can avoid allocating additional memory for the model fields.

Access to index information

PDAF_get_memberid

This routine can be called from the model during an ensemble integration. It provides access to the number (id) of the ensemble member that is currently integrated.

CALL PDAF_get_memberid(memberid)

The only argument is:

  • memberid: In integer providing on output the id the ensemble member

Note:

  • Using PDAF_get_memberid is obviously only useful if more than one ensemble member is integrated by a model task. If there are as many model tasks as ensemble members, memberid is always 1. In this case one can use task_id from the module mod_parallel to distinguish the ensemble members.

PDAF_get_obsmemberid

This routine can be called from the model during the analysis step. It provides access to the number (id) of the ensemble member for which the user-routine for the observation operator is called.

CALL PDAF_get_obsmemberid(memberid)

The only argument is:

  • memberid: In integer providing on output the id the ensemble member

Note:

  • The routine an be useful if the observation operator does not actually operate on the state vector that is provided when obs_op_pdaf is called. Their might be cases in which one likes to read model state information from a file (e.g. if the observation operator performs an averaging over time, while the state vector for the analysis step only contains a single time instance).

Incremental analysis updates

PDAF_incremental

When incremental updating is used, this routine is called during the forecast phase to add a fraction of the analysis increment to an ensemble member.

CALL PDAF_incremental(steps, U_dist_stateinc)

The arguments are:

  • steps: The number of time steps over which the increment should be distributed
  • U_dist_stateinc: The name of the user supplied call-back function that performs the actual addition of the state vector increment to the individual model fields.

PDAF_add_increment

When incremental updating is used, the state increment of the analysis step is not directly added to the forecast state. To add the increment to some state vector PDAF_add_increment is called:

CALL PDAF_add_increment(dim_p, state_p)

The arguments are

  • dim_p: An integer giving the state dimension of state_p
  • state_p: A real array of size dim_p holding the on input the state vector and on output the state vector plus increment.

Note:

  • the routine can be used in the prepoststep routine when the analysis state should be analyzed.

Controlling the MPI communicator for PDAF

PDAF_set_comm_pdaf

This routie allows to sepacify the communicator on which the overall PDAF communication bases.

By default, PDAF bases on MPI_COMM_WORLD, thus all processes in a program. This routine allows to set a different communicator. This can be useful if a model is e.g. run with an OI-server so that the world communicator is split into processes for the IO (file) operations and other processes for the actual model run. In this casem, the model would run using a communicator distinct from MPI_COMM_WORLD and PDAF should operate only with this communicator. PDAF_set_comm_pdaf allows the user to specify this communicator for PDAF.

SUBROUTINE PDAF_set_comm_pdaf(comm_pdaf)

  INTEGER,INTENT(in) :: comm_pdaf    ! MPI communicator for PDAF

Notes:

  • comm_pdaf has to be used consistently in init_parallel_pdaf where the commuicators for PDAF are prepared on the user side.
  • The size of comm_pdaf has to be large enough so that the ensemble run can be performed
  • PDAF_set_comm_pdaf has to be called before calling PDAF_init

Debugging

PDAF_set_debug_flag

This routine was introduced with PDAF C2.1.

the routine allows to activate debugging output. See the documention on PDAF debuggging for more information.

Additional (advanced) functionality

PDAF_reset_forget

This routine was introduced with PDAF V2.0

The routine allows a user to reset the forgetting factor manually during the assimilation process.

For the local ensemble Kalman filters the forgetting factor can be set either globally of differently for each local analysis domain. For the LNETF and the global filters only a global setting of the forgeting factor is possible. This routine allows users to set the forgetting factor case-specific. In addition, the implementation of adaptive choices for the forgetting factor (beyond what is implemented in PDAF) are possible.

The interface is the following:

  SUBROUTINE PDAF_reset_forget(forget)

with the following argument:

  • forget: real, intent(in)
    The new value of the forgetting factor

How to use PDAF_reset_forget

  • Global Filters and LNETF: For global filters PDAF_reset_forget has to be called before the actual analysis step. Within the PDAF context the call can be inserted in prepoststep_pdaf. Alternatively, the call can be inserted in the routine assimilation_pdaf before the call to the PDAF_assimilate_X or PDAF_put_state_X routines
  • Local Filters:
    • global setting: The forgetting factor is reset globally, thus equal for all local analysis dimains, if PDAF_reset_forget is called before the local analysis loop is executed (e.g. in init_obs_pdomi or earlier in the analysis step; see [ImplementAnalysisLocal] for the order of the execution of the call-back routines). Thus, it can also be called in prepoststep_pdaf to get a global seeting as fo rthe global filters.
    • To reset forget specific for each local analysis domain, PDAF_reset_forget should be called during the local analysis loop. In particular the routines init_dim_l or init_dim_obs_l_pdafomi are suited for this.

PDAF_force_analysis

This routine was introduced with PDAF 2.0.

The routine allows a user to enforce the execution of the analysis step at the next call to PDAF_put_state_X or PDAF_assimilate_X.

In particular for PDAF_assimilate_X, the number of time steps is set before the forecast phase is entered. However, one might not know the actual length of the forecast time, e.g. the time when new observation arrive. In this case, one can set for number of time steps to a large value and then check for new observations during the time stepping and call PDAF_force_analysis just before PDAF_assimilate_X is called to enforce that the analysis step is executed.

Warning: This is a routine for advanced functionality of PDAF. Use it carefully, since it can break the assimilation process.

For PDAF_put_state_X, one can use PDAF_force_analysis to overwrite the counting of the ensemble members that PDAF does internally.

  SUBROUTINE PDAF_reset_forget()

without any argument

PDAF_set_memberid

This routine was introduced with PDAF 2.0.

The routine allows a user to reset the ensemble member counter inside PDAF. This routine should only be used in the flexible parallelization mode, thus when PDAF_put_state_X is used. Resetting the member counter, allows e.g. to enforce the execution of the analysis step in the case, when according to the PDAF-internal counter, the ensemble integration is not yet complete. For this, one has to specify the member-ID to be the ensemble size.

To just enforce an analysis step, we recommend to use the routine PDAF_force_analysis. This routine is also compatible with PDAF_assimilate_X.

Warning: This is a routine for advanced functionality of PDAF. Use it carefully, since it can break the assimilation process.

The interface is the following:

  SUBROUTINE PDAF_set_memberid(member)

without the argument:

  • member: integer, intent(in) Member in index to which the PDAF-internal counting is reset.
Last modified 22 months ago Last modified on Feb 22, 2023, 2:04:30 PM