= PDAF_eofcovar = 

This page documents the routine `PDAF_eofcovar` of PDAF, which was introduced with PDAF V1.12.

This routine performs an EOF analysis by singular value decomposition. It is used to prepare a covariance matrix for initializing an ensemble.  For the decomposition a multivariate scaling can be performed to ensure that all fields in the state vectors have unit variance. The EOF decompostion to store a covariance matrix is a convenient way to prepare an ensemble for the data assimilation. One can store the EOFs (singular vectors) and corresponding values in a file. Finally, when one wants to initialize an ensemble from these EOFs for a data assimilation application, one can use the routine [wiki:PDAF_sampleens `PDAF_sampleens`] that generates an ensemble of a chosen size by second-order exact sampling representing the covariance matrix. The ensmeble size can then be chosen freely, but it is limited to the number of EOFs plus one. Thus, it can be useful to stor emore EOFs than one finally might want to use to have the flexibility to cary the ensemble size. 

The routine is typically called in a separate program to prepare and ensembles based on a trajectory of model states.

The interface is the following:
{{{
  SUBROUTINE PDAF_eofcovar(dim, nstates, nfields, dim_fields, offsets, &
       remove_mstate, do_mv, states, stddev, svals, svec, meanstate, status)
}}}
with the following arguments:
 * `dim` : `integer, intent(in)`[[BR]] Dimension of state vector
 * `nstates` : `integer, intent(in)`[[BR]]  Number of state vectors in array `states`
 * `nfields` : `integer, intent(in)`[[BR]] Number of fields in state vector
 * `dim_fields` : `real, intent(in), dimension(nfields)`[[BR]] Size of each field
 * `offsets` : `real, intent(in), dimension(nfields)`[[BR]] Start position of each field in `states`
 * `remove_mstate` : `integer, intent(in)`[[BR]] Set to 1 to compute mean state and subtract it from states
 * `do_mv` : `integer, intent(in)`[[BR]]  Set to 1 to to perform multivariate normalization. `nfields`, `dim_fields`, and `offsets` are only used if `do_mv=1`.
 * `states` : `real, intent(inout), dimension(dim, nstates)`[[BR]] Array of state vectors or state perturbations
 * `stddev` : `real, intent(out), dimension(nfields)`[[BR]] Standard deviations of field variability. Without multivariate scaling (`do_mv=0`), it is `stddev = 1.0`
 * `svals` : `real, intent(out), dimension(nstates)`[[BR]] Singular values divided by sqrt(nstates-1)
 * `svec` : `real, inent(out), dimension(dim, nstates)`[[BR]] Singular vectors
 * `meanstate` : `real, intent(inout), dimension(dim)`[[BR]] Mean state (only used if remove_mstate=1)
 * `status` : `integer, intent(out)`[[BR]] Status flag

To use this routine, one has to initialize the array `states` holding in each column a perturbation vector (state - mean) from a state trajectory. Outputs are the arrays of singular values (`svals`) and left singular vectors (`svec`). The singular values are scaled by '''sqrt(1/(nstates-1))'''. With this '''svec * svals^2^ * svec^T^''' is the covariance matrix. In addition, the standard deviation of the field variations (`stddev`) is an output array, but it is only used if the multivariate normalization is activated (by `do_mv=1`). 

=== Multivariate normalization (do_mv) ===

The multivariate normalization computes for each field in the state vector its standard deviation and then normalizes the field variability for the EOF computation so that it is one. Then the singular value decomposition is computed. Finally the fields are re-scaled to their original values. To use the multivariate normalization one has to define the number of different fields in a state vector (`nfields`), the dimension of all fields (`dim_fields`) and the offset of each field from the start of the state vector (`offsets`). The intention of the multivariate normalization is to ensure that all fields have comparable variability and are hence equally represented by all EOFs. Without the normalization a field with particularly small variability might be essentially absent from the leading EOFs.

The use of the EOF decompositon to generate a covariance matrix with PDAF_eofcovar is exemplified in the Lorenz-96 model example (see testsuite/src/lorenz96/tools/generate_covar.F90).

=== Mean state computation and subtraction (remove_mstate) ===

The option `remove_mstate` allows to tell the routine to compute the mean state of the input states in `states` and to subtract it from this array before computing the EOF decomposition. Finally the mean state is added again to the array `states` to return it unchanged. 

This option allows to let `PDAF_eofcovar` handle the mean-state computation, which is required to compute the EOF decomposition. It is useful for the case that the mean state over all input state has to be computed. If a different mean should be subtracted from the states, the user has to do this before calling `PDAF_eofcovar`. In this case, one set `remove_mstate=0` to avoid a further computation and subtraction of a mean state. This case happens, e.g. if a running mean is subtracted like when the states are snapshots over a model state of the ocean over one year and one wants to remove a running seasonal mean.