Changes between Version 48 and Version 49 of AddFilterAlgorithm


Ignore:
Timestamp:
Nov 7, 2024, 5:29:09 PM (2 weeks ago)
Author:
lnerger
Comment:

--

Legend:

Unmodified
Added
Removed
Modified

  • AddFilterAlgorithm

    v48 v49  
    2121|| `PDAF_options_filters` || interface routine to `PDAF_X_options`.[[BR]] `PDAF_X_options` is an optional routine. Its purpose is to display an overview of available options for the filter algorithm. ||
    2222
     23When `PDAF_init` is called, the DA method is chosen by its ID number. Internally to PDAF, each DA method is identified by a string that is defined in `PDAF_init_filters`. The interface routines have a very simple structure. In general, they select the method-specific routine based on the string identifying the filters. When a DA method is added, a line for the corresponding method-specific routine has to be inserted to each of the interface routines. One can also remove a DA method from PDAF by deleting the corresponding lines form the internal interface routines.
    2324
    2425The routine `PDAF_print_info` only includes the interface to `PDAF_X_memtime`
    2526 * `PDAF_X_memtime` displays information on the run time of the different parts of the assimilation process as well as information on the amount of allocated memory. This functionality is optional.
    2627
    27 
    28 The routines `PDAFomi_assimilate_local` or `PDAFomi_assimilate_global` (or `PDAFomi_put_state_local` or `PDAFomi_put_state_global` for the flexible parallelization variant) are called directly from the model code. These are generic routines which internally call the method-specific routine (PDAF_assimilate_X or PDAF_put_state_X) according to the chosen filter which then calls the actual update routine (PDAF_X_update).
    29 
    30 When `PDAF_init` is called, the filter algorithm is chosen by its ID number. Internally to PDAF, each filter is identified by a string that is defined in `PDAF_init_filters`. The interface routines have a very simple structure. In general, they select the filter-specific routine based on the string identifying the filters. When a filter algorithm is added, a line for the corresponding filter-specific routine has to be inserted to each of the interface routines. One can also remove a filter from PDAF by deleting the corresponding lines form the internal interface routines.
     28The routine `PDAFomi_assimilate_local`, for local ensemble filters and smoothers, or analogously the other available `PDAF_omi_assimilate` routines, like for global ensemble filters and smoothers (`PDAFomi_assimilate_global`) or the different variants of 3D-Var (`PDAFomi_assimilate_*3dvar*`) are called directly from the model code. These are used for the fully-parallel implementation variant. For the flexible-parallelization variant the analogous routines `PDAFomi_put_state_*` are provided. These are generic routines which internally call the method-specific routine (PDAF_assimilate_X or PDAF_put_state_X) according to the chosen filter which then calls the actual update routine (PDAF_X_update). This call structure is explained in Figure 2.
     29
     30[[Image(//pics/analysis_call_structure.png)]]
     31[[BR]]'''Figure 2:''' Internal call structure for the analysis step. The generic routine (here `PDAFomi_assimilate_local` as example) calls the method-specific routine. Here `PDAF_assimilate_X` is the routine that controls the ensemble run in case of the fully-parallel implementation. The routine `PDAF_put_state_X` is used for both the fully-parallel and flexible parallelization variants. For the flexible parallelization one uses `PDAFomi_put_state_local` which directly calls `PDAF_put_state_X`, while for the fully parallel variant this routine is called by `PDAF_assimilate_X`. `PDAF_put_state_X` controls the ensemble for the case that multiple ensmeble states are propagated by a single model task and collects the ensemble from the  ensemble tasks before the assimilate update and distributes the ensemble to the model tasks afterwards. The actual main routine for the DA method is `PDAF_X_update` which is called by `PDAF_put_state_X`. The routines without 'omi' in their name are those will the full interface so that they are usable without OMI.
    3132
    3233== Internal dimensions ==
    … …  
    3536|| `dim_p` || The size of the state vector (with parallelization the size of the local state vector for the current process) ||
    3637|| `dim_ens` || The overall size of the ensemble ||
    37 || `dim_ens_l` || If the ensemble integration is distributed over several ensemble tasks, this variable stores the size of the sub-ensemble handled by the current process. (`dim_ens_l` equals `dim_ens` if no parallelization or if only a single model task is used.)||
    38 || `rank` || The maximum rank of the ensemble covariance matrix. In almost all cases, it is `dim_ens-1`. ||
    39 || `dim_eof` || For mode based filters (currently only SEEK), this is the number of modes used in the state covariance matrix. ||
     38|| `dim_ens_l` || If the ensemble integration is distributed over several ensemble tasks, this variable stores the size of the sub-ensemble handled by the current process. (`dim_ens_l` equals `dim_ens` if no parallelization or if only a single model task is used.) This variable is only used in the ensemble handling of PDAF, but not in the DA update. ||
     39|| `rank` || The maximum rank of the ensemble covariance matrix. In almost all cases, it is `dim_ens-1`. Used in error-subspace filters, (L)ESTKF and (L)SEIK||
    4040
    4141== Internal arrays ==
    4242
    43 Several internal arrays are allocated when PDAF is initialized. These arrays are decraed in `PDAF_mod_filter`. hey are allocated in `PDAF_X_alloc` (see below for details) and remain allocated throughout the assimilation process.
    44 For the processes that compute the analysis (those with `filterpe=.true.`) the following arrays are defined:
     43Several internal arrays are allocated when PDAF is initialized. These arrays are declared in `PDAF_mod_filter`. They are allocated in `PDAF_X_alloc` (see below for details) and remain allocated throughout the assimilation process.
     44For the processes that computes the analysis (those with `filterpe=.true.`) the following arrays are defined:
    4545|| Array || Dimension || Comment ||
    4646||`state`|| `dim_p`   || State vector. Used in all filters. ||
    … …  
    6262 * `PDAF_X_memtime` (optional)
    6363
    64 In addition, the routines
     64The routines are very similar for all DA methods.
     65
     66In addition, one has to implement the routines
    6567 * `PDAF_put_state_X` and `PDAF_assimilate_X`
    66 have to be implemented that are called directly in the model code.
     68
     69As an example we recommend to see e.g. the routines for the ETKF or LETKF.
    6770
    6871We recommend to base on the routines of an existing filter, as most of the routines can be easily adapted to a new filter method.
    … …  
    7679  SUBROUTINE PDAF_X_init(subtype, param_int, dim_pint, param_real, dim_preal, &
    7780                                 ensemblefilter, fixedbasis, verbose, outflag)
    78 }}}
    79 with the following arguments:
    80  * `subtype`: The subtype index of the filter algorithm [integer, input].
    81  * `param_int`: The array of integer parameters [integer(dim_pint), input].
    82  * `dim_pint`: The number of parameters in `param_int` [integer, input].
    83  * `param_real`: The array of real parameters [real(dim_preal), input].
    84  * `dim_preal`: The number of parameters in `param_real` [integer, input]
    85  * `ensemblefilter`: Flag, whether the filter is an ensemble filter or a mode-based filter [logical, output].
    86  * `fixedbasis`: Flag, whether the chosen `subtype` is a filter with fixed ensemble, such that only the ensemble mean is integrated by the model [logical, output].
    87  * `verbose`: Verbosity flag [integer, input]. Valid are the values provided to `PDAF_init`.
    88  * `outflag`: Error flag [integer, output]
    89 
    90 The required operations are to initialize the PDAF-internal parameter variables from the provided values of `subtype`, `param_int`, and `param_real`. In the addition, the logical flags `ensemblefilter` and `fixedbasis` have to be set. The existing implementations also include some screen output about the configuration.
     81
     82  INTEGER, INTENT(inout) :: subtype                ! Sub-type of filter
     83  INTEGER, INTENT(in)    :: dim_pint               ! Number of parameters in param_int
     84  INTEGER, INTENT(inout) :: param_int(dim_pint)    ! The array of integer parameters
     85  INTEGER, INTENT(in)    :: dim_preal              ! Number of parameters in param_real
     86  REAL, INTENT(inout)    :: param_real(dim_preal)  ! The array of real parameters
     87  LOGICAL, INTENT(out)   :: ensemblefilter         ! Flag, whether the filter is an ensemble filter or a mode-based filter
     88  LOGICAL, INTENT(out)   :: fixedbasis             ! Flag, whether the chosen `subtype` uses a fixed ensemble (only the ensemble mean is integrated by the model)
     89  INTEGER, INTENT(in)    :: verbose                ! Control screen output
     90  INTEGER, INTENT(inout) :: outflag                ! Status flag
     91}}}
     92
     93The routine has to perform the following operations:
     94 * initialize the PDAF-internal parameter variables specific for the DA method from the provided values of `subtype`, `param_int`, and `param_real`.
     95 * set the logical flags `ensemblefilter` and `fixedbasis` have to be set.
     96The existing implementations also include some screen output about the configuration.
    9197
    9298Please note:
    9399 * The routine should check, whether the provided value of `subtype` is a valid choice. If this is not the case, the error flag should be set to 2.
    94  * Only parameters from `param_int` and `param_real` up to the value `dim_pint` and `dim_preal` should be considered in the initialization. The array may be bigger, but the user defined which parameters are to be used be setting the values of `dim_pint` and `dim_preal`.
    95  * The error flag `outflag` is initial set to 0.
    96  * The internal parameters are declared and stored in the Fortran module `PDAF_mod_filter`. If a new filter algorithm requires additional parameters, their declaration should be added to the module.
     100 * Only parameters from `param_int` and `param_real` up to the value `dim_pint` and `dim_preal` should be considered in the initialization. The array may be bigger, but the user defines which parameters are to be used be setting the values of `dim_pint` and `dim_preal`.
     101 * The error flag `outflag` is initially set to 0.
     102 * The internal parameters are declared in the Fortran module `PDAF_mod_filter`. If a new filter algorithm requires additional parameters, their declaration should be added to the module. Alternatively, one can introduce a new module the holds the parameters specific to the new DA method.
    97103
    98104
    … …  
    104110{{{
    105111  SUBROUTINE PDAF_X_alloc(subtype, outflag)
    106 }}}
    107 with the following arguments:
    108  * `subtype`: The subtype index of the filter algorithm [integer, input].
    109  * `outflag`: Error flag [integer, input/output]
    110 
    111 All arrays that need to be allocated are declared in the Fortran module `PDAF_mod_filter`. Here, also the dimensions of the arrays are declared. For the allocation of arrays, one has to distinguish between processes that compute the analysis step and those that only participate in the ensemble forecast.
     112
     113  INTEGER, INTENT(in) :: subtype        ! Sub-type of filter
     114  INTEGER, INTENT(out):: outflag        ! Status flag
     115}}}
     116
     117All arrays that need to be allocated are declared in the Fortran module `PDAF_mod_filter`. Here, also the shapes of the arrays are declared. For the allocation of arrays, one has to distinguish between processes that compute the analysis step (`filterpe=.true.`) and those that only participate in the ensemble forecast (`filterpe=.false.`).
    112118
    113119For the processes that compute the analysis (those with `filterpe=.true.`) it is mandatory to allocate the following two arrays:
    … …  
    115121 * `eofV`: This is the ensemble matrix in all ensemble-based filters. For SEEK it is the matrix holding eigenvectors. `eofV` has size (`dim_p`, `dim_ens`).
    116122Depending on the filter algorithm some of the following arrays also need to be allocated:
    117  * `eofU`: This is the eigenvalue matrix '''U''' used in the SEIK and SEEK filters (here, its size is (`rank`,`rank`)). For ETKF, it is the matrix '''A''' of size (`dim_ens`,`dim_ens`). The array only needs to be allocated if the algorithm uses such a matrix. (For EnKF, which does not use this matrix, it is allocated with size (1,1).)
     123 * `eofU`: This is the eigenvalue matrix '''U''' used in the SEIK and SEEK filters and the transform matrix of ESTKF (here, its size is (`rank`,`rank`)). For ETKF, it is the transform matrix '''A''' of size (`dim_ens`,`dim_ens`). The array should always be allocated. For method that don't use such transform matrix, for example the EnKF, it can be allocated with size (1,1).
    118124 * `state_inc`: The increment to the state vector computed in the analysis step. It only needs to be allocated in this routine, if incremental analysis updating is implemented. Otherwise, it is sufficient to allocate and deallocate `state_inc` in the routine performing the analysis step. The size of `state_inc` is `dim_p`.
    119125 * `bias`: If the filter algorithm is implemented with bias correction, the vector `bias` with size `dim_bias_p` is allocated.
    120126
    121 Processes that only participate in the computation of the ensemble forecast, but are not involved in computing the analysis step, operate only on a sub-ensemble. Accordingly, an ensemble array for this sub-ensemble has to be allocated. This is:
     127Processes that only participate in the computation of the ensemble forecast, but are not involved in computing the analysis step (those with `filterpe=.false.`), operate only on a sub-ensemble. Accordingly, an ensemble array for this sub-ensemble has to be allocated. This is:
    122128 * `eofV`: This is the ensemble matrix in all ensemble-based filters. For SEEK it is the matrix holding eigenvectors. For the processes with `filterpe=.false.`, `eofV` has size (`dim_p`, `dim_ens_l`).
    123129
    … …  
    128134The optional routine `PDAF_X_options` displays information on the available options for the filter algorithm.
    129135
     136The rotine has no arguments. Thus the interface is as follows:
     137{{{
     138  SUBROUTINE PDAF_X_options()
     139}}}
     140
     141The following display is recommended:
     142 * Available subtypes (At least '0' for standard implementation)
     143 * Parameters used from the parameter arrays `param_int` and `param_real`.
     144
     145
     146
     147
     148=== `PDAF_X_memtime` ===
     149
     150The optional routine `PDAF_X_memtime` displays information about allocated memory and the execution time of different parts of the filter algorithm.
     151
    130152The interface is as follows:
    131153{{{
    132   SUBROUTINE PDAF_X_options()
    133 }}}
    134 It has no arguments!
    135 
    136 The following display is recommended:
    137  * Available subtypes (At least '0' for standard implementation; '5' for offline mode)
    138  * Parameters used from the parameter arrays `param_int` and `param_real`.
    139 
    140 
    141 
    142 
    143 === `PDAF_X_memtime` ===
    144 
    145 The optional routine `PDAF_X_memtime` displays information about allocated memory and the execution time of different parts of the filter algorithm.
    146 
    147 The interface is as follows:
    148 {{{
    149154  SUBROUTINE PDAF_X_memtime(printtype)
    150 }}}
    151 with the following argument:
    152  * `printtype`: The type of the output to be done [integer, input]. For the filter algorithms that are included in the PDAF source code package the following choices are implemented:
    153    * 1: Display general timers
    154    * 2: Display allocated memory
    155    * 3: Display detailed timers
    156 
    157 The timing operation are implement using the module `PDAF_timer`, which provides the function `PDAF_timeit`. Memory allocation is computed using `PDAF_memcount`, which is provided by the module `PDAF_memcounting`.
    158 
    159 
    160 
    161 === `PDAFomi_put_state_X` / `PDAFomi_assimilate_X` ===
    162 
    163 One of these routines is directly inserted into the model code, if the online mode of PDAF is used. The text on the [ImplementationofAnalysisStep implementation of the analysis step] in the Implementation Guide explains the interface for the algorithms that are included in the PDAF package. Apart from the usual integer status flag, the interface contains the names of the user-supplied routines that are required for the analysis step. Usually, the minimum set of routines are:
    164  * A routine to write model fields back into the ensemble state array (`U_collect_state`).
    165  * A routine that determines the size of the observation vector (`U_init_dim_obs`).
    166  * A routine that contains the implementation of the observation operator (`U_obs_op`).
    167  * The pre- and post-step routine in which the forecast and analysis ensembles can be analyzed (`U_prepoststep`).
    168 Additional routines are possible depending on the requirements of the filter algorithm. When one implements a new filter, one should check, if the filter is compatible with the existing local or global filters. In this case, support for the new filter can be added into one of the existing routines (PDAFomi_assimilate_global/PDAFomi_assimilate_local or likewise PDAFomi_put_state_global/PDAFomi_put_state_local).
    169 For working without OMI one can check if one of the existing routines can be reused. This will facilitate the switching between different filter methods.
    170 
    171 With regard to the operations, `PDAFomi_put_state_X` and `PDAF_put_state_X` prepare for the actual analysis step, that is called inside these routines as a separate routine. The required operations are:
     155
     156  INTEGER, INTENT(in) :: printtype    ! Type of screen output: 
     157                                      ! (1) general timings
     158                                      ! (3) timing information for call-back routines and PDAF-internal operations
     159                                      ! (4) second-level timing information
     160                                      ! (5) very detailed third-level timing information
     161                                      ! (10) process-local allocated memory
     162                                      ! (11) globally allocated memory
     163
     164}}}
     165
     166The timing operation are implemented using the module `PDAF_timer`, which provides the function `PDAF_timeit`. Memory allocation is computed using `PDAF_memcount`, which is provided by the module `PDAF_memcounting`.
     167
     168
     169
     170=== `PDAF_assimilate_X` / `PDAF_put_state_X` ===
     171
     172These routines are called by `PDAFomi_assimilate_*` or `PDAFomi_put_state_*`. They are directly inserted into the model code, if the online mode of PDAF is used. The text on the [ImplementationofAnalysisStep implementation of the analysis step] in the Implementation Guide explains the interface for the algorithms that are included in the PDAF package. As described before `PDAF_assimilate_X` called `PDAF_put_state_X` which then calls the actual DA method. Here `PDAF_assimilate_X` mainly passes its arguments over to `PDAF_put_state_X`.
     173
     174Apart from the usual integer status flag, the interface of the routines contains the names of the user-supplied routines that are required for the analysis step. Usually, the minimum set of routines are:
     175 * `U_collect_state`: The routine that writes model fields into the state vector, i.e. a single column of the ensemble state array
     176 * `U_init_dim_obs`: The routine that determines the size of the observation vector
     177 * `U_obs_op`: The routine that contains the implementation of the observation operator
     178 * `U_init_obs`: The routine that provdes the vector of observations
     179 * `U_prepoststep`: The pre- and post-step routine in which the forecast and analysis ensembles can be analyzed or modified.
     180`PDAF_assimilate_X` uses in addition
     181 * `U_distribute_state`: The routine that fills the field arrays of a model from the state vector, i.e. a single column of the ensemble state array
     182 * `U_next_observation`: The routine that specified the number of time steps until the next DA analysis update.
     183
     184Further routines can be added and depend on the requirements of the filter algorithm. For example the (L)ETKF and (L)ESTKF uses a routine `U_prodRinvA` which has to multiply some temporary matrix of the filter method with the inverse observation error covariance matrix.
     185
     186When one plans to implement a new filter we recommend to check whether the filter is compatible with the existing local or global filters, i.e. using the same call-back routines. In this case, support for the new filter can be added into one of the existing routines (PDAFomi_assimilate_global/PDAFomi_assimilate_local or likewise PDAFomi_put_state_global/PDAFomi_put_state_local).
     187
     188With regard to the operations, the routines `PDAF_put_state_X` prepare for the actual analysis step, that is called inside these routines as a separate routine `PDAF_X_update`. The operations implemented in `PDAF_put_state_X` are:
    172189 * Write model fields back into the ensemble array (by calling `U_collect_state`)
    173190 * Increment the counter for the integrated ensemble members (named `counter` and provided by the module `PDAF_mod_filter`.
    … …  
    178195  * Reset the control variables for the ensemble forecast (`initevol=1`, `member=1`, `step=step_obs+1`).
    179196
    180 In general, the put_state routines of all ensemble-based filters are quite general structures. For the implementation of a new filter one should be able to base on an existing routine, e.g. that of for the ETKF. Then, one has to adapt the interface for the required user-supplied routines of the new filter. In addition, the call of the routine `PDAF_X_update` holding the analysis step has to be revised (name of the routine, required user-supplied routines).
    181 
    182 The routine `PDAF_assimilate_X` is mainy an interface routine to `PDAF_put_state_X`. It counts the time steps and calls `PDAF_put_state_X` when the forecast phase is complete.
     197In general, the PDAF_put_state routines of all ensemble-based filters have the same structure. For the implementation of a new filter we recommend to base on an existing routine, e.g. that of for the ETKF. Thus, one copies the routine to a new name. Then, one adapts the interface for the required user-supplied routines of the new DA method. In addition, the call of the routine `PDAF_X_update` holding the DA analysis update step has to be revised (name of the routine, required user-supplied routines).
     198
     199The routine `PDAF_assimilate_X` is mainy an interface routine to `PDAF_put_state_X`. It counts the time steps and calls `PDAF_put_state_X` when the forecast phase is complete. Thus, apart form adapting the name (`X`) and the interface specifying the call-back routines, there should be no need for changes.
     200
     201== Analysis update step `PDAF_X_update` ==
     202
     203The actual DA analysis update is computed in the routine `PDAF_X_update`. The routine is provided with the names of the call-back routines, the arrays as well of the relevant dimensions as describe above.
     204
     205The structure of the operations in `PDAF_X_update` can be designed freely when implementing a new DA method. The existing methods follow the recommended structure in which `PDAF_x_update` only performes preparations for the actual analysis update and then calls `PDAF_X_analysis` for the computation of the actual update. Here, the operations are different for global and local filters.
     206