Version 25 (modified by 6 weeks ago) (diff)  ,

PDAFOMI: The file callback_obs_pdafomi.F90
PDAFOMI Guide
 Overview
 callback_obs_pdafomi.F90
 Observation Modules
 Observation operators
 Checking error status
 Debugging functionality
 Implementing the analysis step with OMI
 Porting an existing implemention to OMI
 Using domainlimited observations
Contents of this page
The file callback_obs_pdafomi.F90
contains all those routines that are directly called by PDAF as callback 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 passthrough routines. Thus, one typically includes the observationspecific 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:
 Include the observationspecific routine
init_dim_obs_TYPE
and the variableassim_TYPE
from the observationmodule with 'use'  Declare a dimension variable
dim_obs_TYPE
and initialize it to 0  Add a call to the observationspecific routine init_dim_obs_TYPE with the condition
IF (assim_TYPE)
 add
dim_obs_TYPE
to the final sum computingdim_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 ininit_pdaf
, or by reading from the command line or a namelist file.  The obsmodule can have either a specific name for
init_dim_obs_TYPE
or a generic name. If the generic nameinit_dim_obs
is used one can apply a name conversioninit_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 ! PElocal dimension of state INTEGER, INTENT(in) :: dim_obs_p ! Dimension of observed state REAL, INTENT(in) :: state_p(dim_p) ! PElocal model state REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PElocal observed state
The implementation steps are:
 Include the observationspecific routine
obs_op_TYPE
from the observationmodule with 'use'  Add a call to the observationspecific 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 toobs_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 IFstatement with
assim_TYPE
here. The check is done within each obsmodule.
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 domainlocalized filters (like LESTKF, LETKF, LNETF).
The implementation steps are:
 Include the observationspecific routine
init_dim_obs_l_TYPE
from the observationmodule with 'use'  Add a call to the observationspecific 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:
 Include the observationspecific routine
localize_covar_TYPE
from the observationmodule with 'use'  Initialize the array
coords
which holds the coordinates of all elements of the state vector for the current process domain  Add a call to the observationspecific 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 PElocal 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 ininit_pdaf
and declare the array in the modulemod_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 3DVar methods and it is only required in this case (3DVar 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 ! PElocal dimension of state INTEGER, INTENT(in) :: dim_obs_p ! Dimension of observed state REAL, INTENT(in) :: state_p(dim_p) ! PElocal model state REAL, INTENT(out) :: m_state_p(dim_obs_p) ! PElocal 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 toobs_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 IFstatement with
assim_TYPE
here. The check is done within each obsmodule.  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 3DVar methods and it is only required in this case (3DVar 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 ! PElocal dimension of state INTEGER, INTENT(in) :: dim_obs_p ! Dimension of observed state REAL, INTENT(in) :: m_state_p(dim_obs_p) ! PElocal full observed state REAL, INTENT(inout) :: state_p(dim_p) ! PElocal 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 toobs_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 IFstatement with
assim_TYPE
here. The check is done within each obsmodule.
deallocate_obs_pdafomi
The file callback_obs_pdafomi.F90 also contains a routine deallocate_obs_pdafomi. Each obsmodule allocates arrays in the observation type obs_f
and deallocate_obs_pdafomi
is used to deallocate the different observationspecific arrays after the analysis step.
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:
 Include the observationspecific type
thisobs
from the observationmodule with 'use' apply a name conversion likeobs_TYPE => thisobs
 add a call to
PDAFomi_deallocate_obs
giving the observationspecificobs_TYPE
as argument.  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)