PDAF_pf_analysis.F90 Source File


Source Code

! Copyright (c) 2014-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_pf_analysis --- PF analysis with resampling
!
! !INTERFACE:
SUBROUTINE PDAF_pf_analysis(step, dim_p, dim_obs_p, dim_ens, &
     state_p, ens_p, type_resample, noise_type, noise_amp, &
     U_init_dim_obs, U_obs_op, U_init_obs, U_likelihood, &
     screen, flag)

! !DESCRIPTION:
! Analysis step of the PF
!
! Variant for domain decomposed states. 
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2019-05 - Lars Nerger
! 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 PDAF_timer, &
       ONLY: PDAF_timeit
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filtermpi, &
       ONLY: mype
  USE PDAF_mod_filter, &
       ONLY: obs_member
  USE PDAF_mod_filter, &
       ONLY: type_forget, forget, type_winf, limit_winf, debug
  USE PDAFomi, &
       ONLY: omi_n_obstypes => n_obstypes, omi_omit_obs => omit_obs
  USE PDAF_analysis_utils, &
       ONLY: PDAF_omit_obs_omi

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: step         ! Current time step
  INTEGER, INTENT(in) :: dim_p        ! PE-local dimension of model state
  INTEGER, INTENT(out) :: dim_obs_p   ! PE-local dimension of observation vector
  INTEGER, INTENT(in) :: dim_ens      ! Size of ensemble
  REAL, INTENT(inout) :: state_p(dim_p) ! on exit: PE-local forecast mean state
  REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)   ! PE-local state ensemble
  INTEGER, INTENT(in) :: type_resample      ! Type of resampling scheme
  INTEGER, INTENT(in) :: noise_type   ! Type of pertubing noise
  REAL, INTENT(in) :: noise_amp       ! Amplitude of noise
  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_dim_obs, & ! Initialize dimension of observation vector
       U_obs_op, &              ! Observation operator
       U_init_obs, &            ! Initialize observation vector
       U_likelihood             ! Compute observation likelihood for an ensemble member

! !CALLING SEQUENCE:
! Called by: PDAF_pf_update
! Calls: U_init_dim_obs
! Calls: U_obs_op
! Calls: U_init_obs
! Calls: U_likelihood
! Calls: PDAF_timeit
! Calls: PDAF_memcount
! Calls: gemmTYPE (BLAS; dgemm or sgemm dependent on precision)
! Calls: syevTYPE (LAPACK; dsyev or ssyev dependent on precision)
!EOP
       
! *** local variables ***
  INTEGER :: i, col, member           ! counters
  INTEGER, SAVE :: allocflag = 0      ! Flag whether first time allocation is done
  REAL :: effN                        ! Efective sample size
  REAL :: weight                      ! Ensemble weight (likelihood)
  REAL :: total_weight                ! Sum of weights
  INTEGER :: maxblksize, blkupper, blklower  ! Variables for blocked ensemble update
  REAL, ALLOCATABLE :: resid_i(:)     ! PE-local observation residual
  REAL, ALLOCATABLE :: obs_p(:)       ! PE-local observation vector
  REAL, ALLOCATABLE :: Rinvresid(:)   ! R^-1 times residual 
  REAL, ALLOCATABLE :: weights(:)     ! Weight vector
  REAL, ALLOCATABLE :: ens_blk(:,:)   ! Temporary block of state ensemble
  INTEGER, ALLOCATABLE :: IDs(:)      ! Indices for resampled particles


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

  CALL PDAF_timeit(51, 'new')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_pf_analysis -- START'

  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 5x, a)') &
          'PDAF', 'Compute particle filter update'
  END IF


! *********************************
! *** Inflate forecast ensemble ***
! *********************************

  IF (type_forget==0) THEN
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: PDAF_pf_analysis', debug, &
          'Inflate forecast ensemble'
     
     CALL PDAF_timeit(34, 'new') ! Apply forgetting factor
     CALL PDAF_inflate_ens(dim_p, dim_ens, state_p, ens_p, forget)
     CALL PDAF_timeit(34, 'old')
  ENDIF

  CALL PDAF_timeit(51, 'old')


! *********************************
! *** Get observation dimension ***
! *********************************

  IF (debug>0) THEN
     WRITE (*,*) '++ PDAF-debug: PDAF_pf_analysis:', debug, '  dim_p', dim_p
     WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_pf_analysis -- call init_dim_obs'
  END IF


  CALL PDAF_timeit(15, 'new')
  CALL U_init_dim_obs(step, dim_obs_p)
  CALL PDAF_timeit(15, 'old')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_pf_analysis:', debug, '  dim_obs_p', dim_obs_p

  IF (screen > 2) THEN
     WRITE (*, '(a, 5x, a13, 1x, i3, 1x, a, i8)') &
          'PDAF', '--- PE-domain', mype, 'dimension of observation vector', dim_obs_p
  END IF


  ! ***********************************************
  ! *** Compute particle weights (=likelihood)  ***
  ! ***********************************************

  CALL PDAF_timeit(12, 'new')

  ! Allocate weights
  ALLOCATE(weights(dim_ens))   
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens)

  haveobs: IF (dim_obs_p > 0) THEN
     ! *** The weights only exist for domains with observations ***

     ALLOCATE(obs_p(dim_obs_p))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_p)

     ! Allocate tempory arrays for obs-ens_i
     ALLOCATE(resid_i(dim_obs_p))
     ALLOCATE(Rinvresid(dim_obs_p))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', 2*dim_obs_p)

     ! Omit observations if innovation is too large
     ! This step also initializes obs_p
     IF (omi_omit_obs) THEN
       CALL PDAF_omit_obs_omi(dim_p, dim_obs_p, dim_ens, state_p, ens_p, &
            obs_p, U_init_obs, U_obs_op, 1, screen)
     END IF

     ! Get residual as difference of observation and observed state for each ensemble member
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, &
          'PDAF_pf_analysis -- call obs_op and likelihood', dim_ens, 'times'

     CALC_w: DO member = 1, dim_ens

        ! Store member index
        obs_member = member

        CALL PDAF_timeit(44, 'new')
        CALL U_obs_op(step, dim_p, dim_obs_p, ens_p(:, member), resid_i)
        CALL PDAF_timeit(44, 'old')

        IF (member==1 .and. (.not. omi_omit_obs)) THEN
           IF (debug>0) &
                WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_pf_analysis -- call init_obs'

           ! get observation vector (has to be after U_obs_op for OMI)
           CALL PDAF_timeit(50, 'new')
           CALL U_init_obs(step, dim_obs_p, obs_p)
           CALL PDAF_timeit(50, 'old')
        END IF

        CALL PDAF_timeit(51, 'new')
        resid_i = obs_p - resid_i 
        CALL PDAF_timeit(51, 'old')

        IF (debug>0) THEN
           WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_pf_analysis -- member', member
           WRITE (*,*) '++ PDAF-debug PDAF_pf_analysis:', debug, '  innovation d', resid_i
           WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_pf_analysis -- call likelihood'
        end IF

        ! Compute likelihood
        CALL PDAF_timeit(47, 'new')
        CALL U_likelihood(step, dim_obs_p, obs_p, resid_i, weight)
        CALL PDAF_timeit(47, 'old')
        weights(member) = weight

     END DO CALC_w

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_pf_analysis:', debug, '  raw weights', weights

     ! Compute inflation of weights according to N_eff
     IF (type_winf == 1) THEN
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_pf_analysis -- inflate weights'
        CALL PDAF_inflate_weights(screen, dim_ens, limit_winf, weights)
     END IF

     CALL PDAF_timeit(51, 'new')

     ! Normalize weights
     total_weight = 0.0
     DO i = 1, dim_ens
        total_weight = total_weight + weights(i)
     END DO

     IF (total_weight /= 0.0) THEN
        ! Normalize weights
        weights = weights / total_weight

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug PDAF_pf_analysis:', debug, '  normalized weights', weights
     ELSE
        ! weights are zero - reset to uniform weights
        WRITE(*,'(/5x,a/)') 'WARNING: Zero weights - reset to 1/dim_ens'
        weights = 1.0/REAL(dim_ens)
     END IF

     DEALLOCATE(obs_p, resid_i, Rinvresid)

     ! Diagnostic: Compute effective sample size
     CALL PDAF_diag_effsample(dim_ens, weights, effN)
     IF (mype == 0 .AND. screen > 0) &
          WRITE (*, '(a, 5x, a, f10.2)') 'PDAF', '--- Effective sample size ', effN

     CALL PDAF_timeit(51, 'old')

  ELSE
     ! Without observations, all ensemble member have the same weight

     CALL PDAF_timeit(51, 'new')
     weights = 1/dim_ens
     CALL PDAF_timeit(51, 'old')

     ! For OMI we need to call observation operator also for dim_obs_p=0
     ! in order to initialize pointer to observation type
     IF (omi_n_obstypes>0) THEN
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_pf_analysis -- call obs_op', dim_ens, 'times'

        ALLOCATE(resid_i(1))
        obs_member = 1

        ! [Hx_1 ... Hx_N]
        CALL U_obs_op(step, dim_p, dim_obs_p, ens_p(:, 1), resid_i(:))

        DEALLOCATE(resid_i)
     END IF
     
  END IF haveobs

  CALL PDAF_timeit(12, 'old')
  CALL PDAF_timeit(51, 'new')


  ! ****************************************
  ! *** Resample particles               ***
  ! ****************************************

  CALL PDAF_timeit(10, 'new')

  ! Determine sample IDs for resampling

  ALLOCATE(IDs(dim_ens))

  CALL PDAF_timeit(21, 'new')

  CALL PDAF_pf_resampling(type_resample, dim_ens, dim_ens, weights, IDs, screen)

  CALL PDAF_timeit(21, 'old')

  ! Perform the resampling (in blocked form to save memory)

  maxblksize = 200
  IF (mype == 0 .AND. screen > 0) &
       WRITE (*, '(a, 5x, a, i5)') &
       'PDAF', '--- use blocking with size ', maxblksize

  ALLOCATE(ens_blk(maxblksize, dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', maxblksize * dim_ens)

  blocking: DO blklower = 1, dim_p, maxblksize

     blkupper = MIN(blklower + maxblksize - 1, dim_p)

     ! Store old state ensemble
     CALL PDAF_timeit(22, 'new')
     DO col = 1, dim_ens
        ens_blk(1 : blkupper - blklower + 1, col) &
             = ens_p(blklower : blkupper, col)
     END DO

     ! Perform the resampling 
     DO col = 1, dim_ens
        ens_p(blklower : blkupper, col) = ens_blk(1 : blkupper-blklower+1, IDs(col))
     END DO
     CALL PDAF_timeit(22, 'old')

  END DO blocking


  ! *****************************************
  ! *** Perturb particles by adding noise ***
  ! *****************************************

  CALL PDAF_timeit(23, 'new')

  IF (noise_type>0) THEN
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_pf_analysis:', debug, '  add noise to particles'
     CALL PDAF_pf_add_noise(dim_p, dim_ens, state_p, ens_p, noise_type, noise_amp, screen)
  END IF

  CALL PDAF_timeit(23, 'old')


  CALL PDAF_timeit(10, 'old')
  CALL PDAF_timeit(51, 'old')


! *********************************
! *** Inflate analysis ensemble ***
! *********************************

  IF (type_forget==2) THEN
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: PDAF_pf_analysis', debug, &
          'Inflate analysis ensemble'

     CALL PDAF_timeit(34, 'new') ! Apply forgetting factor
     CALL PDAF_inflate_ens(dim_p, dim_ens, state_p, ens_p, forget)
     CALL PDAF_timeit(34, 'old')
  ENDIF


! ********************
! *** Finishing up ***
! ********************
        
  ! Set exit flag
  flag = 0

  DEALLOCATE(ens_blk, IDs)

  IF (allocflag == 0) allocflag = 1

END SUBROUTINE PDAF_pf_analysis