! 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_smoother_lnetf --- Smoother extension for local NETF ! ! !INTERFACE: SUBROUTINE PDAF_smoother_lnetf(domain_p, step, dim_p, dim_l, dim_ens, & dim_lag, Ainv, ens_l, sens_p, cnt_maxlag, & U_g2l_state, U_l2g_state, screen) ! !DESCRIPTION: ! Smoother extension for the ensemble square-root filters (ETKF, ESTKF). ! The routine uses the matrix Ainv computed by the filter analysis ! to perform the smoothing on past ensembles. ! ! Variant for domain decomposed states. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2012-05 - 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 PDAF_timer, & ONLY: PDAF_timeit USE PDAF_memcounting, & ONLY: PDAF_memcount USE PDAF_mod_filtermpi, & ONLY: mype #if defined (_OPENMP) USE omp_lib, & ONLY: omp_get_num_threads, omp_get_thread_num #endif IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: domain_p ! Current local analysis domain INTEGER, INTENT(in) :: step ! Current time step INTEGER, INTENT(in) :: dim_p ! PE-local dimension of model state INTEGER, INTENT(in) :: dim_l ! State dimension on local analysis domain INTEGER, INTENT(in) :: dim_ens ! Size of ensemble INTEGER, INTENT(in) :: dim_lag ! Number of past time instances for smoother REAL, INTENT(in) :: Ainv(dim_ens, dim_ens) ! Weight matrix for ensemble transformation REAL, INTENT(inout) :: ens_l(dim_l, dim_ens) ! local past ensemble (temporary) REAL, INTENT(inout) :: sens_p(dim_p, dim_ens, dim_lag) ! PE-local smoother ensemble INTEGER, INTENT(inout) :: cnt_maxlag ! Count available number of time steps for smoothing INTEGER, INTENT(in) :: screen ! Verbosity flag ! ! External subroutines ! ! (PDAF-internal names, real names are defined in the call to PDAF) EXTERNAL :: U_g2l_state, & ! Get state on local ana. domain from global state U_l2g_state ! Init full state from state on local analysis domain ! !CALLING SEQUENCE: ! Called by: PDAF_etks_update ! Calls: PDAF_timeit ! Calls: PDAF_memcount ! Calls: gemmTYPE (BLAS; dgemm or sgemm dependent on precision) !EOP ! *** local variables *** INTEGER :: member, col, row, lagcol ! Counters INTEGER :: n_lags ! Available number of time instances for smoothing INTEGER :: maxblksize, blkupper, blklower ! Variables for blocked ensemble update INTEGER, SAVE :: allocflag = 0 ! Flag whether first time allocation is done INTEGER, SAVE :: first = 1 ! Flag for very first call to routine INTEGER, SAVE :: domain_save = 1 ! Index of domain from last call to routine REAL :: invdimens ! Inverse of global ensemble size REAL, ALLOCATABLE :: ens_blk(:,:) ! Temporary block of state ensemble REAL, ALLOCATABLE :: W_smooth(:,:) ! Weight matrix for smoothing INTEGER, SAVE :: mythread, nthreads ! Thread variables for OpenMP !$OMP THREADPRIVATE(mythread, nthreads, allocflag, first, domain_save) ! ********************** ! *** INITIALIZATION *** ! ********************** #if defined (_OPENMP) nthreads = omp_get_num_threads() mythread = omp_get_thread_num() #else nthreads = 1 mythread = 0 #endif ! Determine number of time instances for smoothing IF (cnt_maxlag >= dim_lag) THEN ! Already performed enough analysis to smooth over full lag n_lags = dim_lag ELSE ! Not yet enough analysis steps to smoother over full lag n_lags = cnt_maxlag END IF IF ((domain_p <= domain_save) .OR. (first == 1)) THEN IF (mype == 0 .AND. screen > 0 .AND. n_lags > 0 .AND. mythread==0) THEN WRITE (*, '(a, 5x, a, i8)') 'PDAF', 'Perform smoothing up to lag ', n_lags END IF END IF ! ********************************************** ! *** Compute transform matrix for smoothing *** ! *** *** ! *** W_smooth = A_filter + 1_N *** ! ********************************************** havelag: IF (n_lags > 0) THEN invdimens = 1.0 / REAL(dim_ens) ALLOCATE(W_smooth(dim_ens, dim_ens)) IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2) DO col = 1, dim_ens DO row = 1, dim_ens W_smooth(row, col) = Ainv(row, col) END DO END DO ! ********************************************** ! *** Perform smoothing *** ! *** Transform state ensemble at past times *** ! *** a f *** ! *** X = X W_smooth *** ! ********************************************** ! Use block formulation for transformation maxblksize = 200 ALLOCATE(ens_blk(maxblksize, dim_ens)) IF (allocflag == 0) CALL PDAF_memcount(3, 'r', maxblksize * dim_ens) lagcol=1 ! *** Smooth for all available lags *** smoothing: DO lagcol = 1, n_lags ! *** Get local ensemble *** DO member = 1, dim_ens CALL U_g2l_state(step, domain_p, dim_p, sens_p(:, member, lagcol), dim_l, & ens_l(:, member)) END DO ! Use block formulation for transformation blocking: DO blklower = 1, dim_l, maxblksize blkupper = MIN(blklower + maxblksize - 1, dim_l) ! Store former analysis ensemble DO member = 1, dim_ens ens_blk(1 : blkupper-blklower+1, member) & = ens_l(blklower : blkupper, member) END DO ! a f ! Transform ensemble: X = X W_smooth CALL gemmTYPE('n', 'n', blkupper - blklower + 1, dim_ens, dim_ens, & 1.0, ens_blk(1, 1), maxblksize, W_smooth, dim_ens, & 0.0, ens_l(blklower, 1), dim_l) END DO blocking ! *** Initialize global ensemble *** DO member = 1, dim_ens CALL U_l2g_state(step, domain_p, dim_l, ens_l(:, member), dim_p, & sens_p(:, member, lagcol)) END DO END DO smoothing DEALLOCATE(ens_blk, W_smooth) END IF havelag ! ******************** ! *** Finishing up *** ! ******************** ! Increment maxlag counter IF ((domain_p <= domain_save) .OR. (first == 1)) THEN cnt_maxlag = cnt_maxlag + 1 ! Set flag first = 0 END IF domain_save = domain_p ! Set flag for memory counting IF (allocflag == 0) allocflag = 1 END SUBROUTINE PDAF_smoother_lnetf