| 59 | | SUBROUTINE PDAF_assimilate_lenkf(U_collect_state, U_distribute_state, U_init_dim_obs, & |
| 60 | | U_obs_op, U_init_obs, U_prepoststep, U_localize, & |
| 61 | | U_add_obs_err, U_init_obscovar, U_next_observation, status) |
| | 59 | SUBROUTINE PDAF_assimilate_ensrf(U_collect_state, U_distribute_state, & |
| | 60 | U_init_dim_obs_f, U_obs_op_f, U_init_obs_f, U_init_obsvars_f, & |
| | 61 | U_localize_covar_serial, & |
| | 62 | U_prepoststep, U_next_observation, outflag) |
| 64 | | * [#U_collect_statecollect_state_pdaf.F90 U_collect_state]: The name of the user-supplied routine that initializes a state vector from the array holding the ensemble of model states from the model fields. This is basically the inverse operation to `U_distribute_state` used in `PDAF_get_state` as well as here. |
| 65 | | * [#U_distribute_statedistribute_state_pdaf.F90 U_distribute_state]: The name of a user supplied routine that initializes the model fields from the array holding the ensemble of model state vectors. |
| 66 | | * [#U_init_dim_obsinit_dim_obs_pdaf.F90 U_init_dim_obs]: The name of the user-supplied routine that provides the size of observation vector |
| 67 | | * [#U_obs_opobs_op_pdaf.F90 U_obs_op]: The name of the user-supplied routine that acts as the observation operator on some state vector |
| 68 | | * [#U_init_obsinit_obs_pdaf.F90 U_init_obs]: The name of the user-supplied routine that initializes the vector of observations |
| 69 | | * [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep]: The name of the pre/poststep routine as in `PDAF_get_state` |
| 70 | | * [#U_localizelocalize_covar_pdaf.F90 U_localize]: Apply covariance localization to the matrices HP and HPH^T^ |
| 71 | | * [#U_add_obs_erradd_obs_err_pdaf.F90 U_add_obs_err]: The name of the user-supplied routine that adds the observation error covariance matrix to the ensemble covariance matrix projected onto the observation space. |
| 72 | | * [#U_init_obscovarinit_obscovar_pdaf.F90 U_init_obscovar]: The name of the user-supplied routine that initializes the observation error covariance matrix. |
| 73 | | * [#U_next_observationnext_observation_pdaf.F90 U_next_observation]: The name of a user supplied routine that initializes the variables `nsteps`, `timenow`, and `doexit`. The same routine is also used in `PDAF_get_state`. |
| 74 | | * `status`: The integer status flag. It is zero, if `PDAF_assimilate_lenkf` is exited without errors. |
| 75 | | |
| 76 | | == `PDAF_put_state_lenkf` == |
| 77 | | |
| 78 | | When the 'flexible' implementation variant is chosen for the assimilation system, the routine `PDAF_put_state_lenkf` has to be used instead of `PDAF_assimilate_lenkf`. The general aspects of the filter specific routines `PDAF_put_state_*` have been described on the page [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration]. The interface of the routine is identical with that of `PDAF_assimilate_lenkf` with the exception the specification of the user-supplied routines `U_distribute_state` and `U_next_observation` are missing. |
| 79 | | |
| 80 | | The interface when using the LEnKF is the following: |
| 81 | | {{{ |
| 82 | | SUBROUTINE PDAF_put_state_lenkf(U_collect_state, U_init_dim_obs, U_obs_op, & |
| 83 | | U_init_obs, U_prepoststep, U_localize, & |
| 84 | | U_add_obs_err, U_init_obscovar, status) |
| | 65 | * [#U_collect_statecollect_state_pdaf.F90 U_collect_state]: [[BR]]The name of the user-supplied routine that initializes a state vector from the array holding the ensemble of model states from the model fields. This is the inverse operation to `U_distribute_state` used in `PDAF_get_state` as well as here. |
| | 66 | * [#U_distribute_statedistribute_state_pdaf.F90 U_distribute_state]: [[BR]]The name of a user supplied routine that initializes the model fields from the array holding the ensemble of model state vectors. This is the inverse operation to `U_collect_state` |
| | 67 | * [#U_init_dim_obs_finit_dim_obs_f_pdaf.F90 U_init_dim_obs_f]: [[BR]]The name of the user-supplied routine that provides the size of the full observation vector |
| | 68 | * [#U_obs_op_fobs_op_f_pdaf.F90 U_obs_op_f]: [[BR]]The name of the user-supplied routine that acts as the full observation operator on some state vector |
| | 69 | * [#U_init_obs_finit_obs_f_pdaf.F90 U_init_obs_f]: [[BR]]The name of the user-supplied routine that initializes the full vector of observations |
| | 70 | * [#U_init_obsvars_finit_obsvars_pdaf.F90 U_init_obsvars_f]: [[BR]]The name of the user-supplied routine that initializes the vector of observation error variances. |
| | 71 | * [#U_localize_covar_seriallocalize_covar_serial_pdaf.F90 U_localize_covar_serial]: [[BR]]The name of the routine that applies the covariance localization for a single observation |
| | 72 | * [#U_prepoststepprepoststep_ens_pdaf.F90 U_prepoststep]: [[BR]]The name of the pre/poststep routine as in `PDAF_get_state` |
| | 73 | * [#U_next_observationnext_observation_pdaf.F90 U_next_observation]: [[BR]]The name of a user supplied routine that initializes the variables `nsteps`, `timenow`, and `doexit`. The same routine is also used in `PDAF_get_state`. |
| | 74 | * `status_pdaf`: [[BR]]The integer status flag. It is zero, if the routine is exited without errors. |
| | 75 | |
| | 76 | == `PDAF_put_state_ensrf` == |
| | 77 | |
| | 78 | For the offline mode of PDAF, the routine `PDAF_put_state_lenkf` has to be used instead of `PDAF_assimilate_lenkf`. This routine can also be used when the 'flexible' implementation variant is chosen for the assimilation system, The general aspects of the filter specific routines `PDAF_put_state_*` have been described on the page [ModifyModelforEnsembleIntegration Modification of the model code for the ensemble integration]. The interface of the routine is identical with that of `PDAF_assimilate_lenkf` with the exception that the arguments of the user-supplied routines `U_distribute_state` and `U_next_observation` are missing. |
| | 79 | |
| | 80 | The interface is the following: |
| | 81 | {{{ |
| | 82 | SUBROUTINE PDAF_put_state_ensrf(U_collect_state, & |
| | 83 | U_init_dim_obs_f, U_obs_op_f, U_init_obs_f, U_init_obsvars_f, & |
| | 84 | U_localize_covar_serial, & |
| | 85 | U_prepoststep, outflag) |
| 124 | | * It can be useful to not only determine the size of the observation vector at this point. One can also already gather information about the locations of the observations, which will be used later, e.g. to implement the observation operator. An array for the locations can be defined in a module like `mod_assimilation` of the example implementation. |
| 125 | | |
| 126 | | |
| 127 | | === `U_obs_op` (obs_op_pdaf.F90) === |
| 128 | | |
| 129 | | This routine is used by all global filter algorithms (SEEK, SEIK, EnKF, ETKF) and the LEnKF. |
| 130 | | |
| 131 | | The interface for this routine is: |
| 132 | | {{{ |
| 133 | | SUBROUTINE obs_op(step, dim_p, dim_obs_p, state_p, m_state_p) |
| 134 | | |
| 135 | | INTEGER, INTENT(in) :: step ! Currrent time step |
| | 124 | * We recommend to not only determine the size of the observation vector at this point. The routine is a good place to also already gather information about the corresponding indices of the state vector needed later to implement the observation operator. In addition, one can already prepare an array that holds the full observation vector and an array storing the coordinates of the observations. The required arrays can be defined in a module like `mod_assimilation`. The information can be used later in `U_localize_covar_serial`. |
| | 125 | * The routine is similar to `init_dim_obs` used in the global filters. However, if the global filter is used with a domain-decomposed model, it only initializes the size of the observation vector for the local model sub-domain. This is different for the local filters, as the local analysis also requires observational data from neighboring model sub-domains. Nonetheless, one can base on an implemented routine `init_dim_obs` to implement `init_dim_obs_f`. |
| | 126 | |
| | 127 | |
| | 128 | === `U_obs_op_f` (obs_op_f_pdaf.F90) === |
| | 129 | |
| | 130 | This routine is used by the ENSRF and by all filter algorithms with domain-localization (LSEIK, LETKF, LNETF, LKNETF) and is independent of the particular algorithm. |
| | 131 | |
| | 132 | The interface for this routine is: |
| | 133 | {{{ |
| | 134 | SUBROUTINE obs_op_f(step, dim_p, dim_obs_f, state_p, m_state_f) |
| | 135 | |
| | 136 | INTEGER, INTENT(in) :: step ! Current time step |
| 139 | | REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PE-local observed state |
| 140 | | }}} |
| 141 | | |
| 142 | | The routine is called during the analysis step. It has to perform the operation of the observation operator acting on a state vector that is provided as `state_p`. The observed state has to be returned in `m_state_p`. |
| 143 | | |
| 144 | | For a model using domain decomposition, the operation is on the PE-local sub-domain of the model and has to provide the observed sub-state for the PE-local domain. |
| 145 | | |
| 146 | | Hint: |
| 147 | | * If the observation operator involves a global operation, e.g. some global integration, while using domain-decomposition one has to gather the information from the other model domains using MPI communication. |
| 148 | | |
| 149 | | |
| 150 | | === `U_init_obs` (init_obs_pdaf.F90) === |
| 151 | | |
| 152 | | This routine is used by all global filter algorithms (SEEK, SEIK, EnKF, ETKF) and the LEnKF. |
| 153 | | |
| 154 | | The interface for this routine is: |
| 155 | | {{{ |
| 156 | | SUBROUTINE init_obs(step, dim_obs_p, observation_p) |
| 157 | | |
| 158 | | INTEGER, INTENT(in) :: step ! Current time step |
| 159 | | INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of obs. vector |
| 160 | | REAL, INTENT(out) :: observation_p(dim_obs_p) ! PE-local observation vector |
| 161 | | }}} |
| 162 | | |
| 163 | | The routine is called during the analysis step. |
| 164 | | It has to provide the vector of observations in `observation_p` for the current time step. |
| 165 | | |
| 166 | | For a model using domain decomposition, the vector of observations that exist on the model sub-domain for the calling process has to be initialized. |
| | 140 | REAL, INTENT(out) :: m_state_f(dim_obs_f) ! Full observed state |
| | 141 | }}} |
| | 142 | |
| | 143 | The routine is called during the analysis step, before the loop over the single observations is entered. It has to perform the operation of the observation operator acting on a state vector, which is provided as `state_p`. The observed state has to be returned in `m_state_f`. It is the observed state corresponding to the 'full' observation vector. |
| | 144 | |
| | 145 | Hint: |
| | 146 | * The routine is similar to `init_dim_obs` used for the global filters. However, with a domain-decomposed model `m_state_f` will need to contain parts of the state vector from neighboring model sub-domains. Thus, one needs to collect this information which resides in the memory of other processes. PDAF provides the routine [wiki:PDAF_gather_obs_f PDAF_gather_obs_f] for this task. The example implementation in `tutorial/classical/online_2D_parallelmodel` shows the use of `PDAF_gather_obs_f`. |
| | 147 | |
| | 148 | |
| | 149 | === `U_init_obs_f` (init_obs_f_pdaf.F90) === |
| | 150 | |
| | 151 | This routine is used by the ENSRF and by all filter algorithms with domain-localization (LSEIK, LETKF, LNETF, LKNETF) and is independent of the particular algorithm. |
| | 152 | |
| | 153 | The interface for this routine is: |
| | 154 | {{{ |
| | 155 | SUBROUTINE init_obs_f(step, dim_obs_f, observation_f) |
| | 156 | |
| | 157 | INTEGER, INTENT(in) :: step ! Current time step |
| | 158 | INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector |
| | 159 | REAL, INTENT(out) :: observation_f(dim_obs_f) ! Full observation vector |
| | 160 | }}} |
| | 161 | |
| | 162 | The routine is called during the analysis step before the loop over the single observations is entered. It has to provide the full vector of observations in `observation_f` for the current time step. |
| | 163 | |
| | 164 | Hints: |
| | 165 | * As for the other 'full' routines: While the global counterpart of this routine (`init_obs`) has to initialize the observation vector only for the local model sub-domain, the 'full' routine has to include observations that spatially belong to neighboring model sub-domains. As an easy choice one can simply initialize a vector of all globally available observations. |
| | 166 | |
| | 167 | |
| | 168 | === `U_init_obsvars_f` (init_obsvars_f_pdaf.F90) === |
| | 169 | |
| | 170 | This routine is only used by the ENSRF/EAKF. |
| | 171 | |
| | 172 | The interface for this routine is: |
| | 173 | {{{ |
| | 174 | SUBROUTINE init_obsvars_f(step, dim_obs_f, var_f) |
| | 175 | |
| | 176 | INTEGER, INTENT(in) :: step ! Current time step |
| | 177 | INTEGER, INTENT(in) :: dim_obs_f ! Dimension of full observation vector |
| | 178 | REAL, INTENT(out) :: var_f(dim_obs_f) ! vector of observation error variances |
| | 179 | }}} |
| | 180 | |
| | 181 | The routine is called during the analysis step before the loop over the single observations is entered. It has to provide in `var_f` a vector of observation error variances corresponding to the full vector of observations. |
| | 182 | |
| | 183 | === `U_localize_covar_serial` (localize_covar_serial_pdaf.F90) === |
| | 184 | |
| | 185 | This routine is only used by the ENSRF/EAKF. |
| | 186 | |
| | 187 | The interface for this routine is: |
| | 188 | {{{ |
| | 189 | SUBROUTINE U_localize_covar_serial(iobs, dim_p, dim_obs, HP_p, HXY_p) |
| | 190 | |
| | 191 | INTEGER, INTENT(in) :: iobs !< Index of the assimilated single observation |
| | 192 | INTEGER, INTENT(in) :: dim_p !< Process-local state dimension |
| | 193 | INTEGER, INTENT(in) :: dim_obs_f !< Number of full observations |
| | 194 | REAL, INTENT(inout) :: HP_p(dim_p) !< Process-local part of matrix HP for observation iobs |
| | 195 | REAL, INTENT(inout) :: HXY_p(dim_obs_F) !< Process-local part of matrix HX(HX_full) for full observations}}} |
| | 196 | }}} |
| | 197 | |
| | 198 | The routine is called during the loop over all single observations. The purpose of the routine is to apply covariance localization to the vectors '''Hi P''' and '''Hi PH^T^''' for the assimilation of a single observation (determined by the index `iobs` related to the observation operator '''Hi'''). Here '''Hi PH^T^''' is the vector relating to the observed covariance matrix for the full observation vector. This vector is required for parallelization. |
| | 199 | |
| | 200 | Hints: |
| | 201 | * To compute the localization one can use the routine `PDAF_local_weight` after computing the distance between two elements in the vector '''Hi P''' or '''Hi PH^T'''. |
| 206 | | |
| 207 | | === `U_localize` (localize_covar_pdaf.F90) === |
| 208 | | |
| 209 | | This routine is only used for the LEnKF. |
| 210 | | |
| 211 | | The interface for this routine is: |
| 212 | | {{{ |
| 213 | | SUBROUTINE U_localize(dim_p, dim_obs, HP, HPH) |
| 214 | | |
| 215 | | INTEGER, INTENT(in) :: dim_p ! PE-local state dimension |
| 216 | | INTEGER, INTENT(in) :: dim_obs ! Dimension of global observation vector |
| 217 | | REAL, INTENT(inout) :: HP(dim_obs, dim_p) ! Matrix HP |
| 218 | | REAL, INTENT(inout) :: HPH(dim_obs, dim_obs) ! Matrix HPH^T^ |
| 219 | | }}} |
| 220 | | |
| 221 | | The routine is called during the analysis step and has to apply the element-wise Schur product for the covariance localization of the two matrices'''HP''' and '''HPH^T^''', which are provided as input/output arguments. |
| 222 | | |
| 223 | | Notes: |
| 224 | | * In case of a parallelization with domain decomposition, `HP` contains only the columns of the matrix that resides on the model sub-domain of the calling process. The number of rows is that of the global number of observations |
| 225 | | |
| 226 | | Hints: |
| 227 | | * To compute the localization one can use the routine `PDAF_local_weight` after computing the distance between two elements in the matrix '''HP''' or '''HPH^T^'''. |
| 228 | | |
| 229 | | |
| 230 | | === `U_add_obs_err` (add_obs_err_pdaf.F90) === |
| 231 | | |
| 232 | | This routine is only used for the EnKF and LEnKF. |
| 233 | | |
| 234 | | The interface for this routine is: |
| 235 | | {{{ |
| 236 | | SUBROUTINE add_obs_err(step, dim_obs, C) |
| 237 | | |
| 238 | | INTEGER, INTENT(in) :: step ! Current time step |
| 239 | | INTEGER, INTENT(in) :: dim_obs ! Dimension of obs. vector |
| 240 | | REAL, INTENT(inout) :: C(dim_obs, dim_obs) ! Matrix to that the observation |
| 241 | | ! error covariance matrix is added |
| 242 | | }}} |
| 243 | | |
| 244 | | The routine is called during the analysis step. During the analysis step of the LEnKF, the projection of the ensemble covariance onto the observation space is computed. This matrix is provided to the routine as `C_p`. The routine has to add the observation error covariance matrix to `C_p`. |
| 245 | | |
| 246 | | The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization. |
| 247 | | |
| 248 | | Hints: |
| 249 | | * The routine does not require that the observation error covariance matrix is added as a full matrix. If the matrix is diagonal, only the diagonal elements have to be added. |
| 250 | | |
| 251 | | |
| 252 | | |
| 253 | | === `U_init_obscovar` (init_obscovar_pdaf.F90) === |
| 254 | | |
| 255 | | This routine is only used for the EnKF and LEnKF. |
| 256 | | |
| 257 | | The interface for this routine is: |
| 258 | | {{{ |
| 259 | | SUBROUTINE init_obscovar(step, dim_obs, dim_obs_p, covar, m_state_p, & |
| 260 | | isdiag) |
| 261 | | |
| 262 | | INTEGER, INTENT(in) :: step ! Current time step |
| 263 | | INTEGER, INTENT(in) :: dim_obs ! Dimension of observation vector |
| 264 | | INTEGER, INTENT(in) :: dim_obs_p ! PE-local dimension of observation vector |
| 265 | | REAL, INTENT(out) :: covar(dim_obs, dim_obs) ! Observation error covariance matrix |
| 266 | | REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PE-local observation vector |
| 267 | | LOGICAL, INTENT(out) :: isdiag ! Whether the observation error covar. matrix is diagonal |
| 268 | | }}} |
| 269 | | |
| 270 | | The routine is called during the analysis step and is required for the generation of an ensemble of observations. It has to initialize the global observation error covariance matrix `covar`. In addition, the flag `isdiag` has to be initialized to provide the information, whether the observation error covariance matrix is diagonal. |
| 271 | | |
| 272 | | The operation is for the global observation space. Thus, it is independent of whether the filter is executed with or without parallelization. |
| 273 | | |
| 274 | | Hints: |
| 275 | | * The local observation vector `m_state_p` is provided to the routine for the case that the observation errors are relative to the value of the observation. |
| 276 | | |
| 277 | | |
