= PDAF_diag_compute_moments = 

This page documents the routine `PDAF_diag_compute_moments` of PDAF, which was introduced with PDAF V3.0.

This routine computes the mean, the unbiased variance, the unbiased skewness, and the unbiased excess kurtosis from an ensemble. All moments are stored in a single array. `kmax` controls the highest order that should be computed. With parallized models, the moments are computed for the process-local state.

The routine can be called in the pre/poststep routine of PDAF both before and after the analysis step to compute the ensemble statistics.

The interface is the following: 
{{{
SUBROUTINE PDAF_diag_compute_moments(dim_p, dim_ens, ens, kmax, moments, bias) 

  INTEGER, INTENT(IN) :: dim_p              ! Local size of the state vector
  INTEGER, INTENT(IN) :: dim_ens            ! Number of ensemble members/samples
  REAL, INTENT(IN) :: ens(dim_p,dim_ens)    ! Ensemble matrix
  INTEGER, INTENT(IN) :: kmax               ! Highest order of statistical moment (1-4)
  REAL, INTENT(OUT) :: moments(dim_p, kmax) ! The statistical moments up to order `kmax` with
                                            !  - column 1: mean
                                            !  - column 2: variance
                                            !  - column 3: skewness
                                            !  - column 4: excess kurtosis
  INTEGER, OPTIONAL, INTENT(IN) :: bias     !< 0 to compute bias corrected values (default)
}}}


**Note:**
 * The routine does not perform any MPI operations and can also be used if PDAF was not initialized by calling `PDAF_init`.