PDAF_sampleens.F90 Source File


Source Code

! Copyright (c) 2004-2024 Lars Nerger, lars.nerger@awi.de
!
! This routine is free software: you can redistribute it and/or modify
! it under the terms of the GNU Lesser General Public License
! as published by the Free Software Foundation, either version
! 3 of the License, or (at your option) any later version.
!
! PDAF is distributed in the hope that it will be useful,
! but WITHOUT ANY WARRANTY; without even the implied warranty of
! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
! GNU Lesser General Public License for more details.
!
! You should have received a copy of the GNU Lesser General Public
! License along with this software.  If not, see <http://www.gnu.org/licenses/>.
!
!$Id$
!BOP
!
! !ROUTINE: PDAF_SampleEns --- Sample an ensemble from EOF modes
!
! !INTERFACE:
SUBROUTINE PDAF_SampleEns(dim, dim_ens, modes, svals, state, &
     ens, verbose, flag)

! !DESCRIPTION:
! This routine generates an ensemble of model states from a provided
! mean state and EOF modes (singular vectors of a peturbation matrix)
! and singular values. The resulting ensemble is the 2nd order-exact
! sample covariance matrix and mean state. 
! the ensemble state vectors are computed as
!   $ens_i = state + sqrt(dim_ens-1) modes (\Omega C)^T$
! where $C$ holds in its diagonal the singular values ($svals$). $\Omega$
! is an orthogonal transformation matrix that preserves the mean state.
! The generated ensemble fulfills the condition for the state error
! covariance matrix
!   $P = 1/(sqrt(dim_ens-1)  \sum_{i=1}^{dim\_ens} (ens_i - state)(ens_i - state)^T$
!
! !REVISION HISTORY:
! 2014-05 - Lars Nerger - Initial code for SANGOMA based on PDAF example code.
! 2016-05 - Lars Nerger - Back-porting to PDAF
!
! !USES:
! Include definitions for real type of different precision
! (Defines BLAS/LAPACK routines and MPI_REALTYPE)
#include "typedefs.h"

  USE PDAF_mod_filter, &
       ONLY: Nm1vsN, debug

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: dim                   ! Size of state vector
  INTEGER, INTENT(in) :: dim_ens               ! Size of ensemble
  REAL, INTENT(inout) :: modes(dim, dim_ens-1) ! Array of EOF modes
  REAL, INTENT(in)    :: svals(dim_ens-1)      ! Vector of singular values
  REAL, INTENT(inout) :: state(dim)            ! PE-local model state
  REAL, INTENT(out)   :: ens(dim, dim_ens)     ! State ensemble
  INTEGER, INTENT(in) :: verbose               ! Verbosity flag
  INTEGER, INTENT(inout) :: flag               ! Status flag
!EOP

! *** local variables ***
  INTEGER :: row, col                 ! counters
  REAL, ALLOCATABLE :: Omega(:,:)     ! Transformation matrix Omega
  REAL :: fac                         ! Square-root of dim_ens-1


! **********************
! *** INITIALIZATION ***
! **********************

  IF (verbose>0) THEN
     WRITE (*,'(a, 5x,a)') 'PDAF', '******************************************************'
     WRITE (*,'(a, 5x,a)') 'PDAF', '*                  PDAF_SampleEns                    *'
     WRITE (*,'(a, 5x,a)') 'PDAF', '*                                                    *'
     WRITE (*,'(a, 5x,a)') 'PDAF', '*  Sample an ensemble with 2nd-order exact sampling  *'
     WRITE (*,'(a, 5x,a)') 'PDAF', '*    based on the Sangoma tool Sangoma_SampleEns.    *'
     WRITE (*,'(a, 5x,a)') 'PDAF', '******************************************************'

     ! *** Generate full ensemble on filter-PE 0 ***
     WRITE (*, '(/a, 5x, a)') 'PDAF', 'Sample state ensemble from covariance matrix'
     WRITE (*, '(a, 5x, a)') 'PDAF', 'given as EOF vectors and singular values'
     WRITE (*, '(a, 5x, a, i5)') 'PDAF', '--- Ensemble size:  ', dim_ens
     WRITE (*, '(a, 5x, a, i5)') 'PDAF', '--- number of EOFs: ', dim_ens-1
  END IF

  IF (debug>0) THEN
     WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_sampleens -- START'
     WRITE (*,*) '++ PDAF-debug PDAF_samplens:', debug, '  dim', dim
     WRITE (*,*) '++ PDAF-debug PDAF_samplens:', debug, '  modes(1,:)', modes(1,:)
     WRITE (*,*) '++ PDAF-debug PDAF_samplens:', debug, '  svals', svals(:)
     WRITE (*,*) '++ PDAF-debug PDAF_samplens:', debug, &
          '  Note: If REAL values appear incorrect, please check if you provide them with the correct precision'
  END IF

  ! allocate memory for temporary fields
  ALLOCATE(Omega(dim_ens, dim_ens-1))


! ********************************************************
! *** Generate ensemble by transformation of EOF modes ***
! ********************************************************

  ! *** Generate uniform orthogonal matrix OMEGA ***
  CALL PDAF_seik_Omega(dim_ens-1, Omega, 1, verbose)

  ! ***      Generate ensemble of states                  ***
  ! *** ens_i = state + sqrt(dim_ens-1) modes (Omega C)^T ***

  ! A = Omega C
  DO col = 1, dim_ens-1
     DO row = 1, dim_ens
        Omega(row, col) = Omega(row, col) * svals(col)
     END DO
  END DO
      
 ! ens = state + sqrt(dim_ens-1) modes A^T
  DO col = 1, dim_ens
     ens(:, col) = state(:)
  END DO

  IF (Nm1vsN == 1) THEN
     fac = SQRT(REAL(dim_ens-1))
  ELSE
     fac = SQRT(REAL(dim_ens))
  END IF
  CALL gemmTYPE('n', 't', dim, dim_ens, dim_ens-1, &
       fac, modes, dim, Omega, dim_ens, &
       1.0, ens, dim)


! ****************
! *** clean up ***
! ****************

  DEALLOCATE(Omega)

  flag = 0

  IF (debug>0) THEN
     WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_sampleens -- END'
  END IF

END SUBROUTINE PDAF_SampleEns