Changes between Initial Version and Version 1 of OMI_Callback_obs_pdafomi_PDAF23

Timestamp:

Jun 3, 2025, 12:56:16 PM (4 days ago)

Author:

lnerger

Comment:

--

OMI_Callback_obs_pdafomi_PDAF23

@@ +1 @@
+= PDAF-OMI: The file callback_obs_pdafomi.F90 =
+
+{{{
+#!html
+<div class="wiki-toc">
+<h4>PDAF-OMI Guide</h4>
+<ol><li><a href="PDAF_OMI_Overview">Overview</a></li>
+<li>callback_obs_pdafomi.F90</li>
+<li><a href="OMI_observation_modules">Observation Modules</a></li>
+<li><a href="OMI_observation_operators">Observation operators</a></li>
+<li><a href="OMI_error_checking">Checking error status</a></li>
+<li><a href="OMI_debugging">Debugging functionality</a></li>
+<li><a href="OMI_ImplementationofAnalysisStep">Implementing the analysis step with OMI</a></li>
+<ol>
+<li> <a href="ImplementFilterAnalysisOverview"> General overview for ensemble filters</a></li>
+<ol>
+<li><a href="ImplementAnalysisGlobal">Implementation for Global Filters</a></li>
+<li><a href="ImplementAnalysisLocal">Implementation for Local Filters</a></li>
+<li><a href="ImplementAnalysislenkfOmi">Implementation for LEnKF</a></li>
+</ol>
+<li> <a href="Implement3DVarAnalysisOverview"> General overview for 3D-Var methods</a></li>
+<ol>
+<li><a href="ImplementAnalysis_3DVar">Implementation for 3D-Var</a></li>
+<li><a href="ImplementAnalysis_3DEnVar">Implementation for 3D Ensemble Var</a></li>
+<li><a href="ImplementAnalysis_Hyb3DVar">Implementation for Hybrid 3D-Var</a></li>
+</ol>
+</ol>
+<li><a href="OMI_nondiagonal_observation_error_covariance_matrices">Using nondiagonal R-matrices</a></li>
+<li><a href="Porting_to_OMI">Porting an existing implemention to OMI</a></li>
+<li><a href="PDAFomi_additional_functionality">Additional OMI Functionality</a></li>
+</ol>
+</div>
+}}}
+
+[[PageOutline(2-3,Contents of this page)]]
+
+The file `callback_obs_pdafomi.F90` contains all those routines that are directly called by PDAF as call-back routines. In the example codes we use all routines with the suffix _pdafomi to distinguish them from routine from existing implementation where the suffix is typically _pdaf. 
+
+The file `templates/omi/callback_obs_pdafomi.F90` provides a template for implementing the routines. A compact example can be found in `tutorial/online_2D_serialmodel/`.
+
+The routines are mainly pass-through routines. Thus, one typically includes the observation-specific routine with ‘use’ and then calls this routine. However, there is small additional functionality in the different routines which has to be handled when implementing the routine or adding an observation type.
+
+In the descriptions below we use 'TYPE' as a generic label for an observation type. When implemeting an observation type, one replace this by an actual name. For example for sea surface temperature (sst) observations, one could replace TYPE by 'sst' and write e.g. `init_dim_obs_sst`.
+
+
+== init_dim_obs_pdafomi ==
+
+The routine is called at the beginning of each analysis step.  For PDAF, it has to initialize the size `dim_obs_p` of the observation vector according to the current time step. Apart from this routine will initialize overall observation information.  In this routine one just calls `init_dim_obs_TYPE` for each observation type.
+
+The '''interface''' for this routine is:
+{{{
+SUBROUTINE init_dim_obs_pdafomi(step, dim_obs_p)
+
+  INTEGER, INTENT(in)  :: step       ! Current time step
+  INTEGER, INTENT(out) :: dim_obs_p  ! Dimension of observation vector
+}}}
+
+In this routine we declare a dimension variable `dim_obs_TYPE` for each observation type. This is initialized by the corresponding routine `init_dim_obs_TYPE`. The sum of these individual dimensions yields the total number of observations, which is returned to PDAF. 
+
+
+The '''implementation steps''' are:
+1. Include the observation-specific routine `init_dim_obs_TYPE` and the variable `assim_TYPE` from the observation-module with 'use'
+1. Declare a dimension variable `dim_obs_TYPE` and initialize it to 0
+1. Add a call to the observation-specific routine init_dim_obs_TYPE with the condition `IF (assim_TYPE)`
+1. add `dim_obs_TYPE` to the final sum computing `dim_obs`.
+
+As an '''example''', in `tutorial/online_2D_serialmodel/` we have three observations, named A, B, and C. In `init_dim_obs_pdafomi` we find the lines
+{{{
+SUBROUTINE init_dim_obs_pdafomi(step, dim_obs)
+
+  USE obs_A_pdafomi, ONLY: assim_A, init_dim_obs_A
+  USE obs_B_pdafomi, ONLY: assim_B, init_dim_obs_B
+  USE obs_C_pdafomi, ONLY: assim_C, init_dim_obs_C
+
+  ... argument declarations omitted ...
+
+  INTEGER :: dim_obs_A ! Observation dimensions
+  INTEGER :: dim_obs_B ! Observation dimensions
+  INTEGER :: dim_obs_C ! Observation dimensions
+
+  dim_obs_A = 0
+  dim_obs_B = 0
+  dim_obs_C = 0
+
+  IF (assim_A) CALL init_dim_obs_A(step, dim_obs_A)
+  IF (assim_B) CALL init_dim_obs_B(step, dim_obs_B)
+  IF (assim_C) CALL init_dim_obs_C(step, dim_obs_C)
+
+  dim_obs = dim_obs_A + dim_obs_B + dim_obs_C
+}}}
+
+
+'''Notes:'''
+ - The variable `assim_TYPE` indicates whether a particular observation type is assimilated. It is usually set in `init_pdaf`, or by reading from the command line or a namelist file.
+ - The obs-module can have either a specific name for `init_dim_obs_TYPE` or a generic name. If the generic name `init_dim_obs` is used one can apply a name conversion `init_dim_obs_TYPE => init_dim_obs`. 
+
+
+
+== obs_op_pdafomi ==
+
+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`. In this routine one just calls `obs_op_TYPE` for each observation type.
+
+The '''interface''' for this routine is:
+{{{
+SUBROUTINE obs_op_pdafomi(step, dim_p, dim_obs_p, state_p, m_state_p)
+
+  INTEGER, INTENT(in) :: step               ! Currrent time step
+  INTEGER, INTENT(in) :: dim_p              ! PE-local dimension of state
+  INTEGER, INTENT(in) :: dim_obs_p          ! Dimension of observed state
+  REAL, INTENT(in)    :: state_p(dim_p)     ! PE-local model state
+  REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PE-local observed state
+}}}
+
+The '''implementation steps''' are:
+1. Include the observation-specific routine `obs_op_TYPE` from the observation-module with 'use'
+1. Add a call to the observation-specific routine `obs_op_TYPEz
+
+As an '''example''', in `tutorial/online_2D_serialmodel/` we have
+{{{
+SUBROUTINE obs_op_pdafomi(step, dim_p, dim_obs, state_p, ostate)
+
+  USE obs_A_pdafomi, ONLY: obs_op_A
+  USE obs_B_pdafomi, ONLY: obs_op_B
+  USE obs_C_pdafomi, ONLY: obs_op_C
+
+  ... argument declarations omitted ...
+
+  CALL obs_op_A(dim_p, dim_obs, state_p, ostate)
+  CALL obs_op_B(dim_p, dim_obs, state_p, ostate)
+  CALL obs_op_C(dim_p, dim_obs, state_p, ostate)
+}}}
+
+'''Notes:'''
+ - The arguments in the calls to `obs_op_TYPE` are input arguments to `obs_op_pdafomi`. They are just passed on.
+ - The order of the calls to `obs_op_TYPE` is not relevant because the setup of the overall full observation vector is defined by the order of the calls in init_dim_obs_pdafomi. Anyway, it's good practice to keep the order of the calls consistent.
+ - We don't need an IF-statement with `assim_TYPE` here. The check is done within each obs-module.
+
+
+== init_dim_obs_l_pdafomi ==
+
+In this routine one just calls `init_dim_obs_l_TYPE` for each observation type. The routine is only required for domain-localized filters (like LESTKF, LETKF, LNETF).
+
+The '''implementation steps''' are:
+1. Include the observation-specific routine `init_dim_obs_l_TYPE` from the observation-module with 'use'
+1. Add a call to the observation-specific routine init_dim_obs_l_TYPE
+
+As an '''example''', in `tutorial/online_2D_serialmodel/` we have
+{{{
+SUBROUTINE init_dim_obs_l_pdafomi(domain_p, step, dim_obs, dim_obs_l)
+
+  USE obs_A_pdafomi, ONLY: init_dim_obs_l_A
+  USE obs_B_pdafomi, ONLY: init_dim_obs_l_B
+  USE obs_C_pdafomi, ONLY: init_dim_obs_l_C
+
+  ... argument declarations omitted ...
+
+  CALL init_dim_obs_l_A(domain_p, step, dim_obs, dim_obs_l)
+  CALL init_dim_obs_l_B(domain_p, step, dim_obs, dim_obs_l)
+  CALL init_dim_obs_l_C(domain_p, step, dim_obs, dim_obs_l)
+}}}
+
+'''Notes:'''
+ - The order of the calls to `init_dim_obs_l_TYPE` defines the order in which the observations are stored in the local observation vector. The calling order does not need to be the same as in the other routines, but it's good practive to keep the order of the calls consistent.
+
+== localize_covar_pdafomi ==
+
+In this routine one calls `localize_covar_TYPE` for each observation type. The routine is only required for the localized EnKF and performs covariance localization.
+
+The '''implementation steps''' are:
+1. Include the observation-specific routine `localize_covar_TYPE` from the observation-module with 'use'
+1. Initialize the array `coords` which holds the coordinates of all elements of the state vector for the current process domain
+1. Add a call to the observation-specific routine localize_covar_TYPE
+
+As an '''example''', in `tutorial/online_2D_serialmodel/` we have
+{{{
+SUBROUTINE localize_covar_pdafomi(dim_p, dim_obs, HP_p, HPH)
+
+  USE obs_A_pdafomi, ONLY: localize_covar_A
+  USE obs_B_pdafomi, ONLY: localize_covar_B
+  USE obs_C_pdafomi, ONLY: localize_covar_C
+
+  ... argument declarations omitted ...
+
+  REAL, ALLOCATABLE :: coords_p(:,:) ! Coordinates of PE-local state vector entries
+
+  ALLOCATE(coords_p(2, dim_p))
+  coords_p = ...                     ! Initialize coords_p
+
+  CALL localize_covar_A(dim_p, dim_obs, HP_p, HPH, coords_p)
+  CALL localize_covar_B(dim_p, dim_obs, HP_p, HPH, coords_p)
+  CALL localize_covar_C(dim_p, dim_obs, HP_p, HPH, coords_p)
+
+  DEALLOCATE(coords_p)
+}}}
+
+'''Notes:'''
+ - Instead of allocating and filling the coordinate array `coords_p` in this routine one could also do it once in `init_pdaf` and declare the array in the module `mod_assimilation`.
+ - The order of the calls to `obs_op_TYPE` is not relevant because the setup of the overall full observation vector is defined by the order of the calls in init_dim_obs_pdafomi. Anyway, it's good practice to keep the order of the calls consistent.
+
+
+
+== obs_op_lin_pdafomi ==
+
+The routine is called during the analysis step of the 3D-Var methods and it is only required in this case (3D-Var methods have been added with PDAF V2.0). It has to perform the operation of the linearized observation operator acting on a state vector that is provided as `state_p`. The observed state has to be returned in `m_state_p`. In this routine one just calls `obs_op_lin_TYPE` for each observation type. 
+
+The '''interface''' for this routine is:
+{{{
+SUBROUTINE obs_op_lin_pdafomi(step, dim_p, dim_obs_p, state_p, m_state_p)
+
+  INTEGER, INTENT(in) :: step               ! Currrent time step
+  INTEGER, INTENT(in) :: dim_p              ! PE-local dimension of state
+  INTEGER, INTENT(in) :: dim_obs_p          ! Dimension of observed state
+  REAL, INTENT(in)    :: state_p(dim_p)     ! PE-local model state
+  REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PE-local observed state
+}}}
+
+
+The '''implementation steps''' are the same as for `obs_op_pdafomi`.
+
+An '''example''' is provided in `tutorial/variational/online_2D_serialmodel/`.
+
+'''Notes:'''
+ - The arguments in the calls to `obs_op_lin_TYPE` are input arguments to `obs_op_pdafomi`. They are just passed on.
+ - The order of the calls to `obs_op_TYPE` is not relevant because the setup of the overall full observation vector is defined by the order of the calls in init_dim_obs_pdafomi. Anyway, it's good practice to keep the order of the calls consistent.
+ - We don't need an IF-statement with `assim_TYPE` here. The check is done within each obs-module.
+ - If the full observation operator for the observation type is linear by itself, one can just call `obs_op_TYPE`. Thus, no additional routine inside the observation module is required in this case.
+
+
+== obs_op_adj_pdafomi ==
+
+The routine is called during the analysis step of the 3D-Var methods and it is only required in this case (3D-Var methods have been added with PDAF V2.0). It has to perform the operation of the adjoint observation operator acting on an observation vector that is provided as `m_state_p`. The resulting state vector has to be returned in `state_p`. In this routine one just calls `obs_op_adj_TYPE` for each observation type. 
+
+The '''interface''' for this routine is:
+{{{
+SUBROUTINE obs_op_adj_pdafomi(step, dim_p, dim_obs_p, state_p, m_state_p)
+
+  INTEGER, INTENT(in) :: step               ! Currrent time step
+  INTEGER, INTENT(in) :: dim_p              ! PE-local dimension of state
+  INTEGER, INTENT(in) :: dim_obs_p          ! Dimension of observed state
+  REAL, INTENT(in)    :: m_state_p(dim_obs_p) ! PE-local full observed state
+  REAL, INTENT(inout) :: state_p(dim_p)     ! PE-local model state
+}}}
+
+
+The '''implementation steps''' are the same as for `obs_op_pdafomi`.
+
+An '''example''' is provided in `tutorial/variational/online_2D_serialmodel/`.
+
+'''Notes:'''
+ - The arguments in the calls to `obs_op_adj_TYPE` are input arguments to `obs_op_adj_pdafomi`. They are just passed on.
+ - The order of the calls to `obs_op_adj_TYPE` is not relevant because the setup of the overall full observation vector is defined by the order of the calls in init_dim_obs_pdafomi. Anyway, it's good practice to keep the order of the calls consistent.
+ - We don't need an IF-statement with `assim_TYPE` here. The check is done within each obs-module.
+
+
+
+
+
+== deallocate_obs_pdafomi == 
+
+The file callback_obs_pdafomi.F90 also contains a routine deallocate_obs_pdafomi. Each obs-module allocates arrays in the observation type `obs_f` and `deallocate_obs_pdafomi` is used to deallocate the different observation-specific arrays after the analysis step. 
+
+|| [[span(style=color: #FF0000, '''Note:''' Calling `deallocate_obs_pdafomi` is only required in PDAF V1.16. It is no longer required to call it in PDAF V2.0 and later)]] ||
+
+The '''implementation steps''' are:
+1. Include the observation-specific type `thisobs` from the observation-module with 'use' apply a name conversion like `obs_TYPE => thisobs`
+1. add a call to `PDAFomi_deallocate_obs` giving the observation-specific `obs_TYPE` as argument.
+1. To perform the deallocation, insert a call to deallocate_obs_pdafomi at the end of the routine `prepoststep_pdaf`.
+
+As an '''example''', in `tutorial/online_2D_serialmodel/` we have
+{{{
+SUBROUTINE deallocate_obs_pdafomi(step)
+
+  USE PDAFomi, ONLY: PDAFomi_deallocate_obs
+  USE obs_A_pdafomi, ONLY: obs_A => thisobs
+  USE obs_B_pdafomi, ONLY: obs_B => thisobs
+  USE obs_C_pdafomi, ONLY: obs_C => thisobs
+
+  ... argument declarations omitted ...
+
+  CALL PDAFomi_deallocate_obs(obs_A)
+  CALL PDAFomi_deallocate_obs(obs_B)
+  CALL PDAFomi_deallocate_obs(obs_C)
+}}}