PDAF_lknetf_analysis_T.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_lknetf_analysis_T --- LKNETF analysis using T-matrix
!
! !INTERFACE:
SUBROUTINE PDAF_lknetf_analysis_T(domain_p, step, dim_l, dim_obs_l, &
     dim_ens, state_l, Ainv_l, ens_l, HX_l, &
     HXbar_l, state_inc_l, rndmat, forget, &
     obs_l, U_prodRinvA_l, U_init_obsvar_l, U_init_n_domains_p, &
     U_likelihood_l, screen, incremental, type_forget, eff_dimens, &
     type_hyb, hyb_g, hyb_k, gamma, &
     skew_mabs, kurt_mabs, flag)

! !DESCRIPTION:
! Synchronous (1-step) analysis update of the LKNETF. The analysis
! forumulation used a matrix T analogous to the LESTKF and LETKF.
!
! The implementation also supports an adaptive forgetting factor.
!
! Variant for domain decomposed states.
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2017-08 - Lars Nerger - Initial code based on LETKF
! 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: type_trans, obs_member, debug
#if defined (_OPENMP)
  USE omp_lib, &
       ONLY: omp_get_num_threads, omp_get_thread_num
#endif

  IMPLICIT NONE

! !ARGUMENTS:
! ! Variable naming scheme:
! !   suffix _p: Denotes a full variable on the PE-local domain
! !   suffix _l: Denotes a local variable on the current analysis domain
  INTEGER, INTENT(in) :: domain_p    ! Current local analysis domain
  INTEGER, INTENT(in) :: step        ! Current time step
  INTEGER, INTENT(in) :: dim_l       ! State dimension on local analysis domain
  INTEGER, INTENT(in) :: dim_obs_l   ! Size of obs. vector on local ana. domain
  INTEGER, INTENT(in) :: dim_ens     ! Size of ensemble 
  REAL, INTENT(inout) :: state_l(dim_l)           ! local forecast state
  REAL, INTENT(out)   :: Ainv_l(dim_ens, dim_ens) ! on entry: uninitialized
                               ! on exit: local weight matrix for ensemble transformation
  REAL, INTENT(inout) :: ens_l(dim_l, dim_ens)    ! Local state ensemble
  REAL, INTENT(in) :: HX_l(dim_obs_l, dim_ens) ! local observed state ens.
  REAL, INTENT(in) :: HXbar_l(dim_obs_l)       ! local observed ens. mean
  REAL, INTENT(in) :: state_inc_l(dim_l)       ! Local state increment
  REAL, INTENT(inout) :: rndmat(dim_ens, dim_ens) ! Global random rotation matrix
  REAL, INTENT(in) :: obs_l(dim_obs_l) ! Local observation vector
  REAL, INTENT(inout) :: forget        ! Forgetting factor
  REAL, INTENT(inout) :: eff_dimens(1) ! Effective ensemble size
  INTEGER, INTENT(in) :: screen        ! Verbosity flag
  INTEGER, INTENT(in) :: incremental   ! Control incremental updating
  INTEGER, INTENT(in) :: type_forget   ! Type of forgetting factor
  INTEGER, INTENT(in) :: type_hyb      ! Type of hybrid weight
  REAL, INTENT(in) :: hyb_g            ! Prescribed hybrid weight for state transformation
  REAL, INTENT(in) :: hyb_k            ! Scale factor kappa (for type_hyb 3 and 4)
  REAL, INTENT(inout) :: gamma(1)      ! Hybrid weight for state transformation
  REAL, INTENT(inout) :: skew_mabs(1)  ! Mean absolute skewness
  REAL, INTENT(inout) :: kurt_mabs(1)  ! Mean absolute kurtosis
  INTEGER, INTENT(inout) :: flag       ! Status flag

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: &
       U_init_obsvar_l, &      ! Initialize local mean observation error variance
       U_init_n_domains_p, &   ! Provide PE-local number of local analysis domains
       U_prodRinvA_l, &        ! Provide product R^-1 A for local analysis domain
       U_likelihood_l          ! Provide likelihood of an ensemble state

! !CALLING SEQUENCE:
! Called by: PDAF_lknetf_update
! Calls: U_prodRinvA_l
! Calls: U_likelihood_l
! Calls: PDAF_timeit
! Calls: PDAF_memcount
! Calls: PDAF_set_forget_local
! Calls: PDAF_etkf_Tright
! Calls: PDAF_etkf_Tleft
! Calls: PDAF_lknetf_set_gamma
! Calls: gemmTYPE (BLAS; dgemm or sgemm dependent on precision)
! Calls: gemvTYPE (BLAS; dgemv or sgemv dependent on precision)
! Calls: syevTYPE (LAPACK; dsyev or ssyev dependent on precision)
!EOP
       
! *** local variables ***
  INTEGER :: i, j, member, col, row  ! Counters
  INTEGER, SAVE :: allocflag = 0     ! Flag whether first time allocation is done
  INTEGER :: syev_info               ! Status flag for SYEV
  INTEGER :: ldwork                  ! Size of work array for SYEV
  INTEGER :: maxblksize, blkupper, blklower  ! Variables for blocked ensemble update
  REAL    :: sqrtNm1                 ! Temporary variable: sqrt(dim_ens-1)
  REAL :: weight                     ! Ensemble weight (likelihood)
  REAL :: fac                        ! Multiplication factor
  INTEGER, SAVE :: lastdomain = -1   ! store domain index
  LOGICAL :: screenout = .TRUE.      ! Whether to print information to stdout
  REAL, ALLOCATABLE :: HZ_l(:,:)     ! Temporary matrices for analysis
  REAL, ALLOCATABLE :: RiHZ_l(:,:)   ! Temporary matrices for analysis
  REAL, ALLOCATABLE :: resid_l(:)    ! local observation residual
  REAL, ALLOCATABLE :: w_etkf_l(:)   ! local RiHZd
  REAL, ALLOCATABLE :: VRiHZd_l(:)   ! Temporary vector for analysis
  REAL, ALLOCATABLE :: tmp_Ainv_l(:,:) ! Temporary storage of Ainv
  REAL, ALLOCATABLE :: Asqrt_l(:,:)  ! Temporary for square-root of U
  REAL, ALLOCATABLE :: ens_blk(:,:)  ! Temporary block of state ensemble
  REAL, ALLOCATABLE :: svals(:)      ! Singular values of Ainv
  REAL, ALLOCATABLE :: work(:)       ! Work array for SYEV
  REAL, ALLOCATABLE :: resid_i(:)    ! Residual
  REAL, ALLOCATABLE :: weights(:)    ! weight vector
  REAL, ALLOCATABLE :: A_tmp(:,:)    ! Temporary matrix for NETF transformation matrix
  REAL, ALLOCATABLE :: Anetf(:,:)    ! Temporary matrix for NETF transformation matrix
  REAL, ALLOCATABLE :: Asqrt(:,:)    ! Temporary matrix for NETF transformation matrix
  REAL :: total_weight               ! Sum of weight

  INTEGER, SAVE :: mythread, nthreads  ! Thread variables for OpenMP

!$OMP THREADPRIVATE(mythread, nthreads, lastdomain, allocflag, screenout)


! *******************
! *** Preparation ***
! *******************

  CALL PDAF_timeit(51, 'new')

#if defined (_OPENMP)
  nthreads = omp_get_num_threads()
  mythread = omp_get_thread_num()
#else
  nthreads = 1
  mythread = 0
#endif

  ! Control screen output
  IF (lastdomain<domain_p .AND. lastdomain>-1) THEN
     screenout = .FALSE.
  ELSE
     screenout = .TRUE.

     ! In case of OpenMP, let only thread 0 write output to the screen
     IF (mythread>0) screenout = .FALSE.

     ! Output, only in case of OpenMP parallelization
#if defined (_OPENMP)
     IF (screenout) THEN
        WRITE (*,'(a, 5x, a, i5, a)') &
             'PDAF', '--- Use OpenMP parallelization with ', nthreads, ' threads'
     END IF
#endif
  END IF

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

  CALL PDAF_timeit(51, 'old')


! ****************************************************
! *** 1. Weight vector for ensemble mean from ETKF ***
! ****************************************************


! ************************
! *** Compute residual ***
! ***   d = y - H x    ***
! ************************

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

  IF (mype == 0 .AND. screen > 0 .AND. screenout) THEN
     WRITE (*, '(a, 5x, a)') 'PDAF', 'Compute ETKF transform matrix'
  END IF

  haveobsB: IF (dim_obs_l > 0) THEN
     ! *** The residual only exists for domains with observations ***

     ALLOCATE(resid_l(dim_obs_l))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_l)

     ! Get residual as difference of observation and observed state
     resid_l = obs_l - HXbar_l

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  innovation d_l', resid_l

  END IF haveobsB

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


! **********************************************
! ***   Compute analyzed matrix Ainv         ***
! ***                                        ***
! ***     -1                 T  -1           ***
! ***    A  = forget I + (HZ)  R   HZ        ***
! **********************************************

  CALL PDAF_timeit(10, 'new')

  haveobsA: IF (dim_obs_l > 0) THEN
     ! *** The contribution of observation matrix ist only ***
     ! *** computed for domains with observations          ***

     CALL PDAF_timeit(30, 'new')

     ALLOCATE(HZ_l(dim_obs_l, dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_l*dim_ens)
     HZ_l = HX_l

     ! *** Set the value of the forgetting factor  ***
     ! *** Inserted here, because HZ_l is required ***
     CALL PDAF_timeit(51, 'new')
     IF (type_forget == 6) THEN
        CALL PDAF_set_forget_local(domain_p, step, dim_obs_l, dim_ens, HZ_l, &
             HXbar_l, resid_l, obs_l, U_init_n_domains_p, U_init_obsvar_l, &
             forget)
     ENDIF

     ! Subtract ensemble mean: HZ = [Hx_1 ... Hx_N] T
     CALL PDAF_etkf_Tright(dim_obs_l, dim_ens, HZ_l)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF HXT_l', HZ_l(:, 1:dim_ens-1)

     CALL PDAF_timeit(51, 'old')
     CALL PDAF_timeit(30, 'old')
     CALL PDAF_timeit(31, 'new')


     ! ***                RiHZ = Rinv HZ                
     ! *** This is implemented as a subroutine thus that
     ! *** Rinv does not need to be allocated explicitly.
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_analysis -- call prodRinvA_l'

     ALLOCATE(RiHZ_l(dim_obs_l, dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_l * dim_ens)

     CALL PDAF_timeit(48, 'new')
     CALL U_prodRinvA_l(domain_p, step, dim_obs_l, dim_ens, obs_l, HZ_l, RiHZ_l)
     CALL PDAF_timeit(48, 'old')

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF R^-1(HXT_l)', RiHZ_l

     CALL PDAF_timeit(51, 'new')

     ! *** Initialize Ainv = (N-1) I ***
     Ainv_l = 0.0
     DO i = 1, dim_ens
        Ainv_l(i, i) = REAL(dim_ens - 1)
     END DO

     ! ***             T        ***
     ! ***  Compute  HZ  RiHZ   ***

     ALLOCATE(tmp_Ainv_l(dim_ens, dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)

     CALL gemmTYPE('t', 'n', dim_ens, dim_ens, dim_obs_l, &
          1.0, HZ_l, dim_obs_l, RiHZ_l, dim_obs_l, &
          0.0, tmp_Ainv_l, dim_ens)

     DEALLOCATE(HZ_l)
     CALL PDAF_timeit(51, 'old')

  ELSE haveobsA
     ! *** For domains with dim_obs_l=0 there is no ***
     ! *** direct observation-contribution to Ainv  ***

     CALL PDAF_timeit(31, 'new')
     CALL PDAF_timeit(51, 'new')

     ! *** Initialize Ainv = (N-1) I ***
     Ainv_l = 0.0
     DO i = 1, dim_ens
        Ainv_l(i, i) = REAL(dim_ens - 1)
     END DO

     ! No observation-contribution to Ainv from this domain
     ALLOCATE(tmp_Ainv_l(dim_ens, dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)

     tmp_Ainv_l = 0.0

     CALL PDAF_timeit(51, 'old')

  END IF haveobsA


  ! *** Complete computation of Ainv  ***
  ! ***   -1          -1    T         ***
  ! ***  A  =        A  + HZ RiHZ     ***

  CALL PDAF_timeit(51, 'new')
  IF (type_forget<5) THEN
     ! Usually the forgetting factor is not applied here
     Ainv_l = Ainv_l + tmp_Ainv_l
  ELSE
     Ainv_l = forget*Ainv_l + tmp_Ainv_l
  END IF
  CALL PDAF_timeit(51, 'new')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF A^-1_l', Ainv_l

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


! ***********************************************
! *** Compute weight for model state update   ***
! ***                                         ***
! ***              T                    f     ***
! ***    w = A RiHZ d  with d = (y - H x )    ***
! ***                                         ***
! ***********************************************

  CALL PDAF_timeit(13, 'new')
  CALL PDAF_timeit(51, 'new')

  ! *** Compute RiHZd = RiHZ^T d ***
  ALLOCATE(w_etkf_l(dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens)

  haveobsC: IF (dim_obs_l > 0) THEN
     ! *** RiHLd_p/=0 only with observations ***

     CALL gemvTYPE('t', dim_obs_l, dim_ens, 1.0, RiHZ_l, &
          dim_obs_l, resid_l, 1, 0.0, w_etkf_l, 1)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF (HXT_l R^-1)^T d_l', w_etkf_l

     DEALLOCATE(RiHZ_l, resid_l)

  ELSE haveobsC

     w_etkf_l = 0.0

  END IF haveobsC


  ! *** Compute weight vector for state analysis:        ***
  ! ***          w = A RiHZd                             ***
  ! *** Use singular value decomposition of Ainv         ***
  ! ***        Ainv = ASB^T                              ***
  ! *** Then: A = A S^(-1) B                             ***
  ! *** The decomposition is also used for the symmetric ***
  ! *** square-root for the ensemble transformation.     ***

  ! *** Invert Ainv using SVD
  ALLOCATE(svals(dim_ens))
  ALLOCATE(work(3 * dim_ens))
  ldwork = 3 * dim_ens
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', 3 * dim_ens)

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, &
       '  Compute eigenvalue decomposition of ETKF A^-1_l'

  ! Compute SVD of Ainv
  CALL syevTYPE('v', 'l', dim_ens, Ainv_l, dim_ens, svals, work, ldwork, syev_info)

  ! *** check if SVD was successful
  IF (syev_info == 0) THEN
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  eigenvalues', svals

     flag = 0
  ELSE
     WRITE (*, '(/5x, a, i10, a/)') &
          'PDAF-ERROR(1): Domain ', domain_p, ' Problem in SVD of inverse of A_etkf !!!'
     flag = 1
  END IF

  ! *** Compute w = U RiHZd stored in RiHZd
  check0: IF (flag == 0) THEN

     ALLOCATE(VRiHZd_l(dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens)

     CALL gemvTYPE('t', dim_ens, dim_ens, 1.0, Ainv_l, &
          dim_ens, w_etkf_l, 1, 0.0, VRiHZd_l, 1)
     
     DO row = 1, dim_ens
        VRiHZd_l(row) = VRiHZd_l(row) / svals(row)
     END DO
  
     CALL gemvTYPE('n', dim_ens, dim_ens, 1.0, Ainv_l, &
          dim_ens, VRiHZd_l, 1, 0.0, w_etkf_l, 1)

     DEALLOCATE(VRiHZd_l)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF A(HXT_l R^-1)^T d_l', w_etkf_l

  END IF check0

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


! **************************************************************
! *** 2. Weight matrix for ensemble transformation from ETKF ***
! **************************************************************

  CALL PDAF_timeit(51, 'new')

  CALL PDAF_timeit(20, 'new')

  ! Part 1: square-root of U
  DO col = 1, dim_ens
     DO row = 1, dim_ens
        tmp_Ainv_l(row, col) = Ainv_l(row, col) / SQRT(svals(col))
     END DO
  END DO

  ALLOCATE(Asqrt_l(dim_ens, dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)

  sqrtNm1 = SQRT(REAL(dim_ens-1))
  CALL gemmTYPE('n', 't', dim_ens, dim_ens, dim_ens, &
       sqrtNm1, tmp_Ainv_l, dim_ens, Ainv_l, dim_ens, &
       0.0, Asqrt_l, dim_ens)

  ! Optional 
  ! Multiply by orthogonal random matrix with eigenvector (1,...,1)^T
  multrnd: IF (type_trans == 0) THEN
     CALL gemmTYPE('n', 'n', dim_ens, dim_ens, dim_ens, &
          1.0, Asqrt_l, dim_ens, rndmat, dim_ens, &
          0.0, tmp_Ainv_l, dim_ens)
  ELSE
     ! Non-random case
     tmp_Ainv_l = Asqrt_l
  END IF multrnd

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  ETKF A^1/2', tmp_Ainv_l

  CALL PDAF_timeit(20, 'old')


! ****************************************************
! *** 3. Weight vector for ensemble mean from NETF ***
! ****************************************************

  ! **********************************************
  ! *** Compute particle weights as likelihood ***
  ! **********************************************

  IF (mype == 0 .AND. screen > 0 .AND. screenout) THEN
     WRITE (*, '(a, 5x, a)') 'PDAF', 'Compute NETF transform matrix'
  END IF

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

  ! Allocate temporal array for obs-ens_i
  ALLOCATE(resid_i(dim_obs_l))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_l)
     
  CALL PDAF_timeit(22, 'new')
  ! Get residual as difference of observation and observed state for 
  ! each ensemble member only on domains where observations are availible

  CALC_w: DO member = 1,dim_ens

     ! Store current member index
     obs_member = member

     ! Calculate local residual  
     resid_i = obs_l - HX_l(:,member)

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_analysis -- member', member
        WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  innovation d_l', resid_i
        WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_analysis -- call likelihood_l'
     end IF

     ! Compute likelihood
     CALL PDAF_timeit(49, 'new')
     CALL U_likelihood_l(domain_p, step, dim_obs_l, obs_l, resid_i, weight)
     CALL PDAF_timeit(49, 'old')
     weights(member) = weight

  END DO CALC_w

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

  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
     weights = weights / total_weight

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

  DEALLOCATE(resid_i)

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


! **************************************************************
! *** 4. Weight matrix for ensemble transformation from NETF ***
! **************************************************************

  check1: IF (flag == 0) THEN

     ! ****************************************
     ! *** Calculate the transform matrix   ***
     ! ***      A= (diag(w)-w*w^t)          ***
     ! *** with the weights w               ***
     ! ****************************************

     CALL PDAF_timeit(23, 'new')

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

     DO j = 1, dim_ens
        DO i = 1, dim_ens
           Anetf(i,j) = -weights(i) * weights(j)
        ENDDO
     ENDDO
     DO i = 1, dim_ens
        Anetf(i,i) = Anetf(i,i) + weights(i)
     END DO

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  NETF A', Anetf

     CALL PDAF_timeit(23, 'old')


     ! ***************************************
     ! *** Calculate square root of matrix ***
     ! ***************************************

     CALL PDAF_timeit(24, 'new')

     ! Compute symmetric square-root by SVD
     ALLOCATE(A_tmp(dim_ens,dim_ens))
     ALLOCATE(Asqrt(dim_ens,dim_ens))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', 2*dim_ens*dim_ens)
     ldwork = 3*dim_ens
     flag = 0

     CALL PDAF_timeit(31, 'new')

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, &
          '  Compute eigenvalue decomposition of A'

     ! EVD
     CALL syevTYPE('v', 'l', dim_ens, Anetf, dim_ens, svals, work, ldwork, syev_info)

     CALL PDAF_timeit(31, 'old')

     IF (syev_info == 0) THEN
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  eigenvalues', svals
     ELSE
        WRITE(*,'(/5x,a,i7/)') 'Problem computing svd of W-ww^T in domain', domain_p
        flag = 1
     ENDIF
  
     CALL PDAF_timeit(32,'new')  

     ! Ensure to only use positive singular values - negative ones are numerical error
     DO i = 1, dim_ens
        IF (svals(i)>0.0) THEN
           svals(i) = SQRT(svals(i))
        ELSE
           svals(i) = 0.0
        END IF
     END DO

     DO i = 1, dim_ens
        DO j = 1, dim_ens
           Asqrt(j,i) = Anetf(j,i) * svals(i)
        END DO
     END DO

     DEALLOCATE(svals, work)

     CALL PDAF_timeit(32,'old')

     CALL PDAF_timeit(33,'new')

     ! Multiply with singular vectors
     CALL gemmTYPE('n', 't', dim_ens, dim_ens, dim_ens, 1.0, &
          Asqrt, dim_ens, Anetf, dim_ens, 0.0, A_tmp, dim_ens)

     CALL PDAF_timeit(33,'old')

     ! Multiply Asqrt by sqrt(m) to get unbiased ensemble 
     fac = SQRT(REAL(dim_ens))

     CALL PDAF_timeit(35,'new')

     ! Multiply Asqrt with random matrix and the factor 
     CALL gemmTYPE('n', 'n', dim_ens, dim_ens, dim_ens, &
          fac, A_tmp, dim_ens, rndmat, dim_ens, &
          0.0, Asqrt, dim_ens)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  NETF sqrt(N) A^1/2', Asqrt

     CALL PDAF_timeit(35,'old')     

  END IF check1

  
! **********************************
! *** 5. Determine hybrid weight ***
! **********************************

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_analysis -- Compute hybrid gamma'

  CALL PDAF_timeit(53,'new')
  CALL PDAF_lknetf_set_gamma(domain_p, dim_obs_l, dim_ens, &
       HX_l, HXbar_l, weights, type_hyb, hyb_g, hyb_k, &
       gamma, eff_dimens, skew_mabs, kurt_mabs, &
       screen, flag)
  CALL PDAF_timeit(53,'old')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  gamma', gamma


! ************************************************
! ***     6. Transform state ensemble          ***
! ***              a   _f   f                  ***
! ***             X  = X + X  W                ***
! *** The weight matrix W is stored in Ainv_l. ***
! ************************************************


  check2: IF (flag == 0) THEN

     ! **********************************************************
     ! *** Compute transform matrix W including hybridization ***
     ! **********************************************************

     CALL PDAF_timeit(20, 'new')

     IF (mype == 0 .AND. screen > 0 .AND. screenout) THEN
        WRITE (*, '(a, 5x, a)') 'PDAF', 'Perform ensemble transformation'
     END IF

     ! Compute hybrid transformation matrix for ensemble perturbations
     tmp_Ainv_l = gamma(1)*tmp_Ainv_l + (1.0-gamma(1))*Asqrt

     ! Total transformation matrix W = sqrt(U) + w (with hybridization)
     DO col = 1, dim_ens
        DO row = 1, dim_ens
           Ainv_l(row, col) = tmp_Ainv_l(row, col) &
                + gamma(1)*w_etkf_l(row) + (1.0-gamma(1))*weights(row)
        END DO
     END DO

     DEALLOCATE(tmp_Ainv_l, Asqrt, A_tmp)
     DEALLOCATE(w_etkf_l, Asqrt_l, weights)
      
     ! Part 4: T W
     CALL PDAF_etkf_Tleft(dim_ens, dim_ens, Ainv_l)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_analysis:', debug, '  Hybrid transform', Ainv_l

     CALL PDAF_timeit(20, 'old')


     ! ***************************************
     ! *** Perform ensemble transformation ***
     ! ***************************************

     ! Use block formulation for transformation
     maxblksize = 200
     IF (mype == 0 .AND. screen > 0 .AND. screenout) &
          WRITE (*, '(a, 5x, a, i5)') &
          'PDAF', '--- use blocking with size ', maxblksize

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

     blocking: DO blklower = 1, dim_l, maxblksize

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

        ! Store forecast ensemble
        CALL PDAF_timeit(21, 'new')
        DO col = 1, dim_ens
           ens_blk(1 : blkupper - blklower + 1, col) &
                = ens_l(blklower : blkupper, col)
        END DO

        ! Store mean forecast in ensemble matrix
        DO col = 1, dim_ens
           ens_l(blklower : blkupper, col) = state_l(blklower : blkupper)
        END DO
        CALL PDAF_timeit(21, 'old')

        !                        a  _f   f
        ! Transform ensemble:   X = X + X  TW
        CALL PDAF_timeit(22, 'new')

        CALL gemmTYPE('n', 'n', blkupper - blklower + 1, dim_ens, dim_ens, &
             1.0, ens_blk, maxblksize, Ainv_l, dim_ens, &
             1.0, ens_l(blklower, 1), dim_l)

        CALL PDAF_timeit(22, 'old')

     END DO blocking
        
     DEALLOCATE(ens_blk)

  END IF check2

  CALL PDAF_timeit(51, 'old')


! ********************
! *** Finishing up ***
! ********************

  IF (allocflag == 0) allocflag = 1

  ! Store domain index
  lastdomain = domain_p

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_analysis -- END'

END SUBROUTINE PDAF_lknetf_analysis_T