PDAF_set_forget_local.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_set_forget_local - Set local adaptive forgetting factor
!
! !INTERFACE:
SUBROUTINE PDAF_set_forget_local(domain, step, dim_obs_l, dim_ens, mens_l, &
     mstate_l, resid_l, obs_l, U_init_n_domains_p, U_init_obsvar_l, &
     forget)

! !DESCRIPTION:
! Dynamically set the global forgetting factor individually for
! each local analysis domain in LSEIK. 
! This is a typical implementation that tries to ensure
! statistical consistency by enforcing the condition\\
! var\_resid = 1/forget var\_ens + var\_obs\\
! where var\_res is the variance of the innovation residual,
! var\_ens is the ensemble-estimated variance, and
! var\_obs is the observation error variance.\\
! This variant is not proven to improve the estimates!
!
! !REVISION HISTORY:
! 2006-09 - Lars Nerger - Initial code
! Later revisions - see svn log
!
! !USES:
  USE PDAF_timer, &
       ONLY: PDAF_timeit
  USE PDAF_mod_filtermpi, &
       ONLY: mype

  IMPLICIT NONE

! !ARGUMENTS:
  INTEGER, INTENT(in) :: domain           ! Current local analysis domain
  INTEGER, INTENT(in) :: step             ! Current time step
  INTEGER, INTENT(in) :: dim_obs_l        ! Dimension of local observation vector
  INTEGER, INTENT(in) :: dim_ens          ! Ensemble size
  REAL, INTENT(in) :: mens_l(dim_obs_l, dim_ens) ! Local observed ensemble
  REAL, INTENT(in) :: mstate_l(dim_obs_l) ! Local observed state estimate
  REAL, INTENT(in) :: resid_l(dim_obs_l)  ! Local observation-forecast residual
  REAL, INTENT(in) :: obs_l(dim_obs_l)    ! Local observation vector
  REAL, INTENT(inout) :: forget           ! Forgetting factor

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: U_init_n_domains_p, & ! Provide number of local analysis domains
       U_init_obsvar_l              ! Initialize local mean obs. error variance

! !CALLING SEQUENCE:
! Called by: PDAF_lseik_analysis
! Calls: U_init_n_domains_p
! Calls: U_init_obsvar_l
!EOP
  
! *** local variables ***
  INTEGER :: i, j                     ! Counters
  REAL :: var_ens, var_resid, var_obs ! Variances
  INTEGER, SAVE :: first = 1          ! Flag for very first call to routine
  INTEGER, SAVE :: domain_save = 1    ! Index of domain from last call to routine
  INTEGER :: n_domains                ! Number of local analysis domains
  REAL :: forget_neg, forget_max, forget_min ! limiting values of forgetting factor


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

  ! Define limiting values of forgetting factor
  ! These are set very arbitrarily for now
  forget_neg = forget
  forget_max = 100.0
  forget_min = 0.01


! ****************************************************
! *** Initialize adaptive local forgetting factors ***
! ****************************************************

  ! Get number of local analysis domains
  CALL PDAF_timeit(42, 'new')
  CALL U_init_n_domains_p(step, n_domains)
  CALL PDAF_timeit(42, 'old')

  IF ((domain <= domain_save) .OR. (first == 1)) THEN
     ! At first call during each forecast phase
     IF (mype == 0) THEN
        WRITE (*, '(a, 5x, a)') &
             'PDAF', '--- Apply dynamically estimated local forgetting factors'
        WRITE (*, '(a, 9x, a, es10.2)') &
             'PDAF', 'Maximum limit for forgetting factor', forget_max
        WRITE (*, '(a, 9x, a, es10.2)') &
             'PDAF', 'Minimum limit for forgetting factor', forget_min
        WRITE (*, '(a, 9x, a, es10.2)') &
             'PDAF', 'Forgetting factor if var(obs)>var(resid)', forget_neg
     ENDIF
          
     ! Set flag
     first = 0
  ENDIF
  domain_save = domain


  ! ************************************************************
  ! *** Compute optimal forgetting factor for current domain ***
  ! ************************************************************

  ! *** Compute mean ensemble variance ***

  ! local
  var_ens = 0.0
  DO i = 1, dim_obs_l
     DO j = 1, dim_ens
        var_ens = var_ens + (mstate_l(i) - mens_l(i, j)) ** 2
     ENDDO
  ENDDO
  var_ens = var_ens / REAL(dim_ens - 1) / REAL(dim_obs_l)


  ! *** Compute mean of observation-minus-forecast residual ***
   
  var_resid = 0.0
  DO i = 1, dim_obs_l
     var_resid = var_resid + resid_l(i) ** 2
  ENDDO
  var_resid = var_resid / REAL(dim_obs_l)


  ! *** Compute mean observation variance ***

  ! Get mean observation error variance
  CALL PDAF_timeit(52, 'new')
  CALL U_init_obsvar_l(domain, step, dim_obs_l, obs_l, var_obs)
  CALL PDAF_timeit(52, 'old')

  ! *** Compute optimal forgetting factor ***
  forget = var_ens / (var_resid - var_obs)

  ! Apply special condition if observation variance is larger than residual variance
  IF (forget < 0.0) forget = forget_neg

  ! Impose upper limit for forgetting factor
  ! - the value for this is quite arbitary
  IF (forget > forget_max) forget = forget_max

  ! Impose lower limit for forgetting factor
  ! - the value for this is quite arbitary
  IF (forget < forget_min) forget = forget_min
   
END SUBROUTINE PDAF_set_forget_local