! 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 - Set adaptive forgetting factor ! ! !INTERFACE: SUBROUTINE PDAF_set_forget(step, filterstr, dim_obs_p, dim_ens, mens_p, & mstate_p, obs_p, U_init_obsvar, forget_in, forget_out) ! !DESCRIPTION: ! Dynamically set the global forgetting factor. ! 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 routine is used in SEIK. It can also be used in LSEIK. ! In this case a forgetting factor for the PE-local domain is ! computed. An alternative for LSEIK is PDAF\_set\_forget\_local, ! which computes a forgetting factor for each local analysis ! domain. The implementation used in both routines is ! experimental and not proven to improve the estimates. ! ! !REVISION HISTORY: ! 2006-09 - 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_mod_filtermpi, & ONLY: mype, npes_filter, MPIerr, COMM_filter, dim_eof_l IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: step ! Current time step CHARACTER(len=10), INTENT(in) :: filterstr ! String defining filter INTEGER, INTENT(in) :: dim_obs_p ! Dimension of observation vector INTEGER, INTENT(in) :: dim_ens ! Ensemble size REAL, INTENT(in) :: mens_p(dim_obs_p, dim_eof_l) ! Observed PE-local ensemble REAL, INTENT(in) :: mstate_p(dim_obs_p) ! Observed PE-local mean state REAL, INTENT(in) :: obs_p(dim_obs_p) ! Observation vector REAL, INTENT(in) :: forget_in ! Prescribed forgetting factor REAL, INTENT(out) :: forget_out ! Adaptively estimated forgetting factor ! ! External subroutine ! ! (PDAF-internal name, real name is defined in the call to PDAF) EXTERNAL :: U_init_obsvar ! Initialize mean obs. error variance ! !CALLING SEQUENCE: ! Called by: PDAF_seik_analysis, PDAF_seik_analysis_newT ! Called by: PDAF_lseik_update ! Calls: U_init_obsvar ! Calls: MPI_allreduce (MPI) !EOP ! *** local variables *** INTEGER :: i, j ! Counters INTEGER :: dim_obs ! Global observation dimension REAL :: var_ens_p, var_ens ! Variance of ensemble REAL :: var_resid_p, var_resid ! Variance of residual REAL :: var_obs ! Variance of observation errors 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_in forget_max = 100.0 forget_min = 0.01 IF (mype == 0) THEN WRITE (*, '(a, 5x, a)') & 'PDAF', '--- Apply global adaptive forgetting factor' 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 ! ****************************************** ! *** Compute adaptive forgetting factor *** ! ****************************************** haveobs: IF (dim_obs_p > 0) THEN ! *** Compute mean ensemble variance *** CALL PDAF_timeit(51, 'new') IF (TRIM(filterstr) /= 'LSEIK' .AND. TRIM(filterstr) /= 'LETKF') THEN ! global IF (npes_filter>1) THEN CALL MPI_allreduce(dim_obs_p, dim_obs, 1, & MPI_INTEGER, MPI_SUM, COMM_filter, MPIerr) ELSE ! This is a work around for working with nullmpi.F90 dim_obs = dim_obs_p END IF ELSE ! For LSEIK use only PE-local observation dimension dim_obs = dim_obs_p ENDIF ! local var_ens_p = 0.0 DO i = 1, dim_obs_p DO j = 1, dim_ens var_ens_p = var_ens_p + (mstate_p(i) - mens_p(i, j)) ** 2 ENDDO ENDDO var_ens_p = var_ens_p / REAL(dim_ens - 1) / REAL(dim_obs) IF (TRIM(filterstr) /= 'LSEIK' .AND. TRIM(filterstr) /= 'LETKF') THEN ! global CALL MPI_allreduce(var_ens_p, var_ens, 1, MPI_REALTYPE, MPI_SUM, & COMM_filter, MPIerr) ELSE ! For LSEIK use only PE-local variance var_ens = var_ens_p ENDIF ! *** Compute mean of innovation *** ! Compute variance var_resid_p = 0.0 DO i = 1, dim_obs_p var_resid_p = var_resid_p + (obs_p(i) - mstate_p(i)) ** 2 ENDDO var_resid_p = var_resid_p / REAL(dim_obs) IF (TRIM(filterstr) /= 'LSEIK' .AND. TRIM(filterstr) /= 'LETKF') THEN ! global CALL MPI_allreduce(var_resid_p, var_resid, 1, & MPI_REALTYPE, MPI_SUM, COMM_filter, MPIerr) ELSE ! For LSEIK use only PE-local variance var_resid = var_resid_p ENDIF CALL PDAF_timeit(51, 'old') ! *** Compute mean observation variance *** ! Get mean observation error variance CALL PDAF_timeit(49, 'new') CALL U_init_obsvar(step, dim_obs_p, obs_p, var_obs) CALL PDAF_timeit(49, 'old') CALL PDAF_timeit(51, 'new') ! *** Compute optimal forgetting factor *** forget_out = var_ens / (var_resid - var_obs) ! Apply special condition if observation variance is larger than residual variance IF (forget_out < 0.0) forget_out = forget_neg ! Impose upper limit for forgetting factor IF (forget_out > forget_max) forget_out = forget_max ! Impose lower limit for forgetting factor IF (forget_out < forget_min) forget_out = forget_min ! ******************** ! *** FINISHING UP *** ! ******************** IF (mype == 0) THEN WRITE (*, '(a, 9x, a, es10.2)') & 'PDAF', '--> Computed forgetting factor', forget_out ENDIF CALL PDAF_timeit(51, 'old') ENDIF haveobs END SUBROUTINE PDAF_set_forget