PDAF_enkf_obs_ensemble.F90 Source File


Source Code

! Copyright (c) 2004-2024 Lars Nerger
!
! This file is part of PDAF.
!
! PDAF 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 PDAF.  If not, see <http://www.gnu.org/licenses/>.
!
!$Id$
!BOP
!
! !ROUTINE: PDAF_enkf_obs_ensemble --- Generate ensemble of observations for EnKF
!
! !INTERFACE:
SUBROUTINE PDAF_enkf_obs_ensemble(step, dim_obs_p, dim_obs, dim_ens, m_ens_p, &
     U_init_obs, U_init_obs_covar, screen, flag)

! !DESCRIPTION:
! This routine generates an ensemble of observations 
! from a mean observation with prescribed error
! covariance matrix for the EnKF94/98.
! For a diagonal matrix the column vectors are 
! directly used. For a non-diagonal matrix the 
! eigen-decomposition is computed.
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2001-01 - Lars Nerger - Initial code
! Later revisions - see svn log
!
! !USES:
! Include definitions for real type of different precision
! (Defines BLAS/LAPACK routines and MPI_REALTYPE)
#include "typedefs.h"

  USE mpi
  USE PDAF_timer, &
       ONLY: PDAF_timeit
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filtermpi, &
       ONLY: mype, npes_filter, MPIerr, COMM_filter
  USE PDAF_mod_filter, &
       ONLY: debug
  USE PDAFomi, &
       ONLY: omi_n_obstypes => n_obstypes, map_obs_id

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: step       ! Current time step
  INTEGER, INTENT(in) :: dim_obs_p  ! Local dimension of current observation
  INTEGER, INTENT(in) :: dim_obs    ! PE-local dimension of observation vector
  INTEGER, INTENT(in) :: dim_ens    ! Size of ensemble
  REAL, INTENT(out)   :: m_ens_p(dim_obs_p,dim_ens) ! PE-local obs. ensemble 
  INTEGER, INTENT(in) :: screen     ! Verbosity flag
  INTEGER, INTENT(inout) :: flag    ! Status flag

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: U_init_obs, & ! Initialize observation vector
       U_init_obs_covar     ! Initialize observation error covariance matrix

! !CALLING SEQUENCE:
! Called by: PDAF_enkf_analysis_rlm
! Called by: PDAF_enkf_analysis_rsm
! Calls: U_init_obs
! Calls: U_init_obs_covar
! Calls: syevTYPE (LAPACK; dsyev or ssyev dependent on precision)
! Calls: larnvTYPE (BLAS; dlarnv or slarnv dependent on precision)
! Calls: MPI_allgather (MPI)
!EOP

! *** local variables ***
  INTEGER :: i, j, member         ! Counters
  REAL :: randval                 ! Value of random number
  REAL, ALLOCATABLE :: m_state_p(:) ! Observation vector
  REAL, ALLOCATABLE :: covar(:, :)  ! Observation covariance matrix
  INTEGER :: syev_info            ! Output flag of eigenproblem routine
  INTEGER, SAVE :: allocflag = 0  ! Flag for first-time allocation
  INTEGER, SAVE :: iseed(4)       ! Seed for random number generator LARNV
  INTEGER, SAVE :: first = 1      ! Flag for setting of random-number seed
  LOGICAL :: isdiag               ! Is the observation error cov. matrix diagonal?
  REAL, ALLOCATABLE :: eigenv(:)  ! Vector of eigenvalues
  REAL, ALLOCATABLE :: workarray(:) ! Workarray for eigenproblem routine
  INTEGER, ALLOCATABLE :: local_dim_obs(:) ! Array of local dimensions
  INTEGER, ALLOCATABLE :: local_dis(:)     ! Array of local displacements
  REAL, ALLOCATABLE :: randvals(:)  ! Vector of random numbers


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

  CALL PDAF_timeit(51, 'new')

  IF (mype == 0 .AND. screen > 0) &
       WRITE (*, '(a, 5x, a)') 'PDAF', '--- Generate ensemble of observations'

  IF (first == 1) THEN
     ! Initialize seed
     iseed(1) = 1
     iseed(2) = 5
     iseed(3) = 7
     iseed(4) = 9
     first = 2
  END IF

  ! allocate memory for temporary fields
  ALLOCATE(eigenv(dim_obs))
  ALLOCATE(workarray(3 * dim_obs))
  ALLOCATE(m_state_p(dim_obs_p))
  ALLOCATE(covar(dim_obs, dim_obs))
  IF (allocflag == 0) THEN
     ! count allocated memory
     CALL PDAF_memcount(3, 'r', 4 * dim_obs + dim_obs_p + dim_obs*dim_obs)
     allocflag = 1
  END IF

  CALL PDAF_timeit(51, 'old')


! *************************************
! *** generate observation ensemble ***
! *************************************

  ! *** get current observation vector ***
  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_enkf_obs_ensemble -- call init_obs'

  CALL PDAF_timeit(50, 'new')
  CALL U_init_obs(step, dim_obs_p, m_state_p)
  CALL PDAF_timeit(50, 'old')

  ! *** Get current observation covariance matrix ***
  ! *** We initialize the global observation error covariance matrix
  ! *** to avoid a parallelization of the possible eigendecomposition.
  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_enkf_obs_ensemble -- call init_obs_covar'

  CALL PDAF_timeit(49, 'new')
  covar = 0.0
  CALL U_init_obs_covar(step, dim_obs, dim_obs_p, covar, m_state_p, &
       isdiag)
  CALL PDAF_timeit(49, 'old')

  CALL PDAF_timeit(51, 'new')

  diagA: IF (.NOT. isdiag) THEN
     ! *** compute Eigendecomposition of covariance matrix ***
     ! *** We use the LAPACK routine SYEV                  ***
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_enkf_obs_ensemble:', debug, &
          '  Compute eigenvalue decomposition of cvoarance matrix R'

     CALL syevTYPE('v', 'l', dim_obs, covar, dim_obs, &
          eigenv, workarray, 3 * dim_obs, syev_info)

     IF (syev_info==0 .AND. debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_enkf_resample:', debug, &
          '  eigenvalues (1:min(dim_obs,20))', eigenv(1:min(dim_obs,20))

  ELSE diagA
     ! *** Do not perform eigendecomposition here, since COVAR is diagonal
     IF (mype == 0 .AND. screen > 0) &
         WRITE (*, '(a, 5x, a)') 'PDAF', '--- use diagonal observation eror cov. matrix'
     syev_info = 0
  END IF diagA

  ! check if eigendecomposition was successful
  ensemble: IF (syev_info /= 0) THEN
     ! Eigendecomposition failed

     WRITE (*, '(/5x, a/)') &
          'PDAF-ERROR(3): Problem in eigendecomposition of observation covariance !!!'
     flag = 3

  ELSE
     ! Eigendecomposition OK, continue with ensemble generation

     diagB: IF (.NOT. isdiag) THEN
        ! rescale eigenvectors (if EVP was computed)
        DO j = 1, dim_obs
           DO i = 1, dim_obs
              covar(i, j) = covar(i, j) * SQRT(eigenv(j))
           END DO
        END DO
     ELSE diagB
        ! Only compute square-roots of variances here
        DO i = 1, dim_obs
           covar(i, i) = SQRT(covar(i, i))
        END DO
     END IF diagB

     ! gather array of observation dimensions
     ALLOCATE(local_dim_obs(npes_filter))
     ALLOCATE(local_dis(npes_filter))

     CALL MPI_allgather(dim_obs_p, 1, MPI_INTEGER, local_dim_obs, 1, &
          MPI_INTEGER, COMM_filter, MPIerr)

     ! Init array of displacements
     local_dis(1) = 0
     DO i = 2, npes_filter
        local_dis(i) = local_dis(i - 1) + local_dim_obs(i - 1)
     END DO

     USE_OMI: IF (omi_n_obstypes == 0) THEN
        ! If not using OMI

        ! generate random states for local domain
        members: DO member = 1, dim_ens
           m_ens_p(:, member) = m_state_p(:)
           eigenvectors: DO j = 1, dim_obs
              CALL larnvTYPE(3, iseed, 1, randval)
              components: DO i = 1, dim_obs_p
                 m_ens_p(i, member) = m_ens_p(i, member) &
                      + covar(i + local_dis(mype + 1), j) * randval
              END DO components
           END DO eigenvectors
        END DO members


     ELSE USE_OMI
        ! For OMI use the matting vector map_obs_id to ensure consistency
        ! if different numbers of processes are used.

        ALLOCATE(randvals(dim_obs))

        ! generate random states for local domain
        membersB: DO member = 1, dim_ens

           m_ens_p(:, member) = m_state_p(:)

           ! Create vector of random numbers
           DO j = 1, dim_obs
              CALL larnvTYPE(3, iseed, 1, randvals(j))
           END DO

           ! Create perturbed observations
           eigenvectorsB: DO j = 1, dim_obs
              componentsB: DO i = 1, dim_obs_p
                 m_ens_p(i, member) = m_ens_p(i, member) &
                      + covar(i + local_dis(mype + 1), j) * randvals(map_obs_id(j))
              END DO componentsB
           END DO eigenvectorsB
        END DO membersB

        DEALLOCATE(randvals)
     END IF USE_OMI

     DEALLOCATE(local_dim_obs, local_dis)

  END IF ensemble

  CALL PDAF_timeit(51, 'old')

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

  DEALLOCATE(m_state_p, covar)
  DEALLOCATE(eigenv, workarray)

END SUBROUTINE PDAF_enkf_obs_ensemble