PDAF_lknetf_step_update.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_step_update --- Control analysis update of the 2-step LKNETF
!
! !INTERFACE:
SUBROUTINE  PDAF_lknetf_step_update(step, dim_p, dim_obs_f, dim_ens, &
     state_p, Uinv, ens_p, state_inc_p, &
     U_init_dim_obs, U_obs_op, U_init_obs, U_init_obs_l, U_prodRinvA_hyb_l, &
     U_init_n_domains_p, U_init_dim_l, U_init_dim_obs_l, U_g2l_state, U_l2g_state, &
     U_g2l_obs, U_init_obsvar, U_init_obsvar_l, U_likelihood_l, U_likelihood_hyb_l, &
     U_prepoststep, screen, subtype, incremental, type_forget, &
     dim_lag, sens_p, cnt_maxlag, flag)

! !DESCRIPTION:
! Routine to control the analysis update of the 2-step LKNETF.
!
! The analysis is performed by first preparing several
! global quantities on the PE-local domain, like the
! observed part of the state ensemble for all local
! analysis domains on the PE-local state domain.
! Then the analysis is done in two steps local analysis 
! domains in the PE-local state domain.  Either the NETF is 
! computed before the LETKF (subtype=0) or LETKF is 
! computed before NETF (subtype=1). 
! In the local analysis loops, the local state and 
! observation dimensions are initialized and the global 
! state ensemble is restricted to the local analysis domain.
! In addition, the routine U\_prepoststep is called prior
! to the analysis and after the resampling outside of
! the loop over the local domains to allow the user
! to access the ensemble information.
!
! !  This is a core routine of PDAF and
!    should not be changed by the user   !
!
! !REVISION HISTORY:
! 2018-01 - Lars Nerger - Initial code based on PDAF_lknetf_update
! 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, PDAF_time_temp
  USE PDAF_memcounting, &
       ONLY: PDAF_memcount
  USE PDAF_mod_filter, &
       ONLY: obs_member, type_trans, filterstr, forget, inloop, &
       member_save, type_hyb, hyb_g, hyb_k, &
       skewness, kurtosis, store_rndmat, debug
  USE PDAFomi, &
       ONLY: omi_omit_obs => omit_obs
  USE PDAF_mod_filtermpi, &
       ONLY: mype, dim_ens_l, npes_filter, COMM_filter, MPIerr
  USE PDAF_analysis_utils, &
       ONLY: PDAF_print_domain_stats, PDAF_init_local_obsstats, PDAF_incr_local_obsstats, &
       PDAF_print_local_obsstats

  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
! !   suffix _f: Denotes a full variable of all observations required for the
! !              analysis loop on the PE-local domain
  INTEGER, INTENT(in) :: step          ! Current time step
  INTEGER, INTENT(in) :: dim_p         ! PE-local dimension of model state
  INTEGER, INTENT(out) :: dim_obs_f    ! PE-local dimension of observation vector
  INTEGER, INTENT(in) :: dim_ens       ! Size of ensemble
  REAL, INTENT(inout) :: state_p(dim_p)         ! PE-local model state
  REAL, INTENT(inout) :: Uinv(dim_ens, dim_ens) ! Inverse of matrix U
  REAL, INTENT(inout) :: ens_p(dim_p, dim_ens)  ! PE-local ensemble matrix
  REAL, INTENT(inout) :: state_inc_p(dim_p)     ! PE-local state analysis increment
  INTEGER, INTENT(in) :: screen        ! Verbosity flag
  INTEGER, INTENT(in) :: subtype       ! Filter subtype
  INTEGER, INTENT(in) :: incremental   ! Control incremental updating
  INTEGER, INTENT(in) :: type_forget   ! Type of forgetting factor
  INTEGER, INTENT(in) :: dim_lag       ! Number of past time instances for smoother
  REAL, INTENT(inout) :: sens_p(dim_p, dim_ens, dim_lag) ! PE-local smoother ensemble
  INTEGER, INTENT(inout) :: cnt_maxlag ! Count number of past time steps for smoothing
  INTEGER, INTENT(inout) :: flag       ! Status flag

! ! External subroutines 
! ! (PDAF-internal names, real names are defined in the call to PDAF)
  EXTERNAL :: U_obs_op, &    ! Observation operator
       U_init_n_domains_p, & ! Provide number of local analysis domains
       U_init_dim_l, &       ! Init state dimension for local ana. domain
       U_init_dim_obs, &     ! Initialize dimension of observation vector
       U_init_dim_obs_l, &   ! Initialize dim. of obs. vector for local ana. domain
       U_init_obs, &         ! Initialize observation vector
       U_init_obs_l, &       ! Init. observation vector on local analysis domain
       U_init_obsvar, &      ! Initialize mean observation error variance
       U_init_obsvar_l, &    ! Initialize local mean observation error variance
       U_g2l_state, &        ! Get state on local ana. domain from global state
       U_l2g_state, &        ! Init full state from state on local analysis domain
       U_g2l_obs, &          ! Restrict full obs. vector to local analysis domain
       U_prodRinvA_hyb_l, &  ! Compute product of R^(-1) with HV with hybrid weight
       U_likelihood_l, &     ! Compute likelihood
       U_likelihood_hyb_l, & ! Compute likelihood with hybrid weight
       U_prepoststep         ! User supplied pre/poststep routine

! !CALLING SEQUENCE:
! Called by: PDAF_put_state_lknetf
! Calls: U_prepoststep
! Calls: U_init_n_domains_p
! Calls: U_init_dim_obs
! Calls: U_obs_op
! Calls: U_init_obs
! Calls: U_init_dim_l
! Calls: U_init_dim_obs_l
! Calls: U_g2l_state
! Calls: U_l2g_state
! Calls: U_g2l_obs
! Calls: PDAF_set_forget
! Calls: PDAF_inflate_ens
! Calls: PDAF_generate_rndmat
! Calls: PDAF_lknetf_compute_gamma
! Calls: PDAF_lknetf_ana_lnetf
! Calls: PDAF_lknetf_ana_letkfT
! Calls: PDAF_timeit
! Calls: PDAF_memcount
! Calls: MPI_reduce
!EOP

! *** local variables ***
  INTEGER :: i, member, row        ! Counters
  INTEGER :: domain_p              ! Counter for local analysis domain
  INTEGER, SAVE :: allocflag = 0   ! Flag whether first time allocation is done
  REAL    :: invdimens             ! Inverse global ensemble size
  INTEGER :: minusStep             ! Time step counter
  INTEGER :: n_domains_p           ! number of PE-local analysis domains
  REAL    :: forget_ana_l          ! forgetting factor supplied to analysis routine
  REAL    :: forget_ana            ! Possibly globally adaptive forgetting factor
  LOGICAL :: storerndmat = .FALSE. ! Store and reuse random rotation matrix
  REAL, ALLOCATABLE :: HX_f(:,:)   ! HX for PE-local ensemble
  REAL, ALLOCATABLE :: HXbar_f(:)  ! PE-local observed mean state
  REAL, ALLOCATABLE :: obs_f(:)    ! PE-local observation vector
  REAL, ALLOCATABLE :: rndmat(:,:) ! random rotation matrix for ensemble trans.
  REAL, SAVE, ALLOCATABLE :: rndmat_save(:,:) ! Stored rndmat
  REAL :: invforget                ! inverse forgetting factor
  ! Variables on local analysis domain
  INTEGER :: dim_l                 ! State dimension on local analysis domain
  INTEGER :: dim_obs_l             ! Observation dimension on local analysis domain
  REAL, ALLOCATABLE :: HX_l(:,:)   ! HX for local analysis domain
  REAL, ALLOCATABLE :: HXbar_l(:)  ! local observed mean state
  REAL, ALLOCATABLE :: obs_l(:)    ! local observation vector
  REAL, ALLOCATABLE :: resid_l(:)  ! local residual
  REAL, ALLOCATABLE :: ens_l(:,:)  ! State ensemble on local analysis domain
  REAL, ALLOCATABLE :: state_l(:)  ! Mean state on local analysis domain
  REAL, ALLOCATABLE :: stateinc_l(:)  ! State increment on local analysis domain
  REAL, ALLOCATABLE :: n_eff(:)    ! Effective sample size for each local domain (in alpha)
  REAL, ALLOCATABLE :: n_eff_b(:)  ! Effective sample size for each local domain (hybrid)
  LOGICAL, ALLOCATABLE :: MASK(:)  ! Mask for effective sample sizes > 0
  REAL :: max_n_eff_l, min_n_eff_l ! PE-local min/max. effective ensemble sizes
  REAL :: max_n_eff, min_n_eff     ! Global min/max. effective ensemble sizes
  REAL :: max_gamma_l, min_gamma_l ! PE-local min/max. hybrid weight
  REAL :: max_gamma, min_gamma     ! Global min/max. hybrid weight
  REAL :: sum_gamma_l, mean_gamma  ! Local alpha sum; global mean alpha
  REAL :: sum_n_eff_l, mean_n_eff  ! Local sum of N_eff; global mean N_eff
  REAL :: max_stats_l(2), min_stats_l(2)   ! PE-local min/max of skewness and kurtosis
  REAL :: max_stats(2), min_stats(2)  ! Global min/max of skewness and kurtosis
  REAL :: sum_stats_l(2)           ! PE-local sum of skewness and kurtosis for averaging
  REAL :: mean_stats(2)            ! Global average skewness and kurtosis
  REAL, ALLOCATABLE :: Uinv_l(:,:) ! thread-local matrix Uinv
  REAL :: state_inc_p_dummy        ! Dummy variable to avoid compiler warning
  REAL, ALLOCATABLE :: gamma(:)    ! Hybrid weight for state update
  INTEGER :: cnt_small_svals       ! Counter for small values
  INTEGER :: n_domains_with_obs_p  ! Domain-local number of local domains with observations
  INTEGER :: n_domains_with_obs    ! Global number of local domains with observations


! ***********************************************************
! *** For fixed error space basis compute ensemble states ***
! ***********************************************************

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

  ! Initialize variable to prevent compiler warning
  state_inc_p_dummy = state_inc_p(1)


! *************************************
! *** Prestep for forecast ensemble ***
! *************************************

  CALL PDAF_timeit(5, 'new')
  minusStep = - step  ! Indicate forecast by negative time step number
  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 5x, a, i7)') 'PDAF', 'Call pre-post routine after forecast; step ', step
  ENDIF
  CALL U_prepoststep(minusStep, dim_p, dim_ens, dim_ens_l, dim_obs_f, &
       state_p, Uinv, ens_p, flag)
  CALL PDAF_timeit(5, 'old')

  IF (mype == 0 .AND. screen > 0) THEN
     IF (screen > 1) THEN
        WRITE (*, '(a, 5x, a, F10.3, 1x, a)') &
             'PDAF ', '--- duration of prestep:', PDAF_time_temp(5), 's'
     END IF
     WRITE (*, '(a, 55a)') 'PDAF Analysis ', ('-', i = 1, 55)
  END IF


! **************************************
! *** Preparation for local analysis ***
! **************************************

#ifndef PDAF_NO_UPDATE
  CALL PDAF_timeit(3, 'new')
  CALL PDAF_timeit(4, 'new')

  IF (debug>0) THEN
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_int(3) -not used-     '
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_int(4) -not used-  '
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_int(5) type_forget', type_forget
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_int(6) type_trans ', type_trans
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_int(7) type_hyb   ', type_hyb

     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_real(1) forget   ', forget
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_real(2) hyb_gamma', hyb_g
     WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update', debug, &
          'Configuration: param_real(3) hyb_kappa', hyb_k
  END IF

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_n_domains'

  ! Query number of analysis domains for the local analysis
  ! in the PE-local domain
  CALL PDAF_timeit(42, 'new')
  CALL U_init_n_domains_p(step, n_domains_p)
  CALL PDAF_timeit(42, 'old')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  n_domains_p', n_domains_p

  ! Initialize effective sample sizes and MASK
  ALLOCATE(n_eff(n_domains_p))
  ALLOCATE(n_eff_b(n_domains_p))
  n_eff = 0.0
  n_eff_b = 0.0
  ALLOCATE(MASK(n_domains_p))

  ! Allocate arrays for skewness and kurtosis
  IF(ALLOCATED(skewness)) DEALLOCATE(skewness)
  IF(ALLOCATED(kurtosis)) DEALLOCATE(kurtosis)
  ALLOCATE(skewness(n_domains_p))
  ALLOCATE(kurtosis(n_domains_p))
  skewness = 0.0
  kurtosis = 0.0

  ! Allocate arrays for hybrid weights
  ALLOCATE(gamma(n_domains_p))
  gamma = 0.0
  
  IF (screen > 0) THEN
     IF (mype == 0) THEN
        IF (subtype == 0) THEN
           WRITE (*, '(a, i7, 3x, a)') &
                'PDAF ', step, 'Assimilating observations - 2-step LKNETF-HNK: NETF before LETKF'
        ELSE IF (subtype == 1) THEN
           WRITE (*, '(a, i7, 3x, a)') &
                'PDAF ', step, 'Assimilating observations - 2-step LKNETF-HKN: LETKF before NETF'
        END IF
     END IF
     IF (screen<3) THEN
        CALL PDAF_print_domain_stats(n_domains_p)
     ELSE
        WRITE (*, '(a, 5x, a, i6, a, i10)') &
             'PDAF', '--- PE-domain:', mype, ' number of analysis domains:', n_domains_p
     END IF
  END IF


! *** Local analysis: initialize global quantities ***

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_dim_obs'

  ! Get observation dimension for all observations required 
  ! for the loop of local analyses on the PE-local domain.
  CALL PDAF_timeit(43, 'new')
  CALL U_init_dim_obs(step, dim_obs_f)
  CALL PDAF_timeit(43, 'old')

  IF (screen > 2) THEN
     WRITE (*, '(a, 5x, a, i6, a, i10)') &
          'PDAF', '--- PE-Domain:', mype, &
          ' dimension of PE-local full obs. vector', dim_obs_f
  END IF

  ! Apply covariance inflation to global ens (only if prior infl is chosen)
  ! returns the full ensemble with unchanged mean, but inflated covariance
  IF (type_forget==0 .OR. type_forget==1) THEN
     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- apply inflation to forecast ensemble'

     CALL PDAF_timeit(14, 'new')
     CALL PDAF_inflate_ens(dim_p, dim_ens, state_p, ens_p, forget)
     CALL PDAF_timeit(14, 'old')
  ENDIF

  ! HX = [Hx_1 ... Hx_(r+1)] for full DIM_OBS_F region on PE-local domain
  ALLOCATE(HX_f(dim_obs_f, dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_f * dim_ens)

  CALL PDAF_timeit(44, 'new')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call obs_op', dim_ens, 'times'

  ENS: DO member = 1,dim_ens
     ! Store member index to make it accessible with PDAF_get_obsmemberid
     obs_member = member

     ! Call observation operator
     CALL U_obs_op(step, dim_p, dim_obs_f, ens_p(:, member), HX_f(:, member))
  END DO ENS

  CALL PDAF_timeit(44, 'old')

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

  ! *** Compute mean forecast state
  state_p = 0.0
  invdimens = 1.0 / REAL(dim_ens)
  DO member = 1, dim_ens
     DO row = 1, dim_p
        state_p(row) = state_p(row) + invdimens * ens_p(row, member)
     END DO
  END DO

  CALL PDAF_timeit(11, 'old')

  ! *** Compute mean state of ensemble on PE-local observation space 
  ALLOCATE(HXbar_f(dim_obs_f))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_f)

  HXbar_f = 0.0
  invdimens = 1.0 / REAL(dim_ens)
  DO member = 1, dim_ens
     DO row = 1, dim_obs_f
        HXbar_f(row) = HXbar_f(row) + invdimens * HX_f(row, member)
     END DO
  END DO
  CALL PDAF_timeit(51, 'old')

  ! Set forgetting factor globally
  forget_ana = forget
  IF (type_forget == 5) THEN
     ALLOCATE(obs_f(dim_obs_f))
     IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_obs_f)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_obs_op'

     ! get observation vector
     CALL PDAF_timeit(50, 'new')
     CALL U_init_obs(step, dim_obs_f, obs_f)
     CALL PDAF_timeit(50, 'old')

     ! Set FORGET
     CALL PDAF_set_forget(step, filterstr, dim_obs_f, dim_ens, HX_f, &
          HXbar_f, obs_f, U_init_obsvar, forget, forget_ana)
     
     DEALLOCATE(obs_f)
  ENDIF


  ! *** Initialize random transformation matrix
  CALL PDAF_timeit(51, 'new')
  CALL PDAF_timeit(33, 'new')
  ALLOCATE(rndmat(dim_ens, dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)

  rnd_store: IF (.NOT. storerndmat .OR. (storerndmat .AND. allocflag == 0)) THEN

     IF (type_trans == 0) THEN
        ! Initialize random matrix
        IF (screen > 0 .AND. mype == 0) &
             WRITE (*, '(a, 6x, a)') 'PDAF', '--- Initialize random transformation'
        CALL PDAF_generate_rndmat(dim_ens, rndmat, 2)

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  rndmat', rndmat
     ELSE
        IF (screen > 0 .AND. mype == 0) &
             WRITE (*, '(a, 6x, a)') 'PDAF', '--- Initialize deterministic transformation'
        rndmat = 0.0
        DO i = 1, dim_ens
           rndmat(i,i) = 1.0
        END DO
     END IF

     IF (store_rndmat) THEN
        IF (allocflag == 0) ALLOCATE(rndmat_save(dim_ens, dim_ens))
        IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)
        rndmat_save = rndmat
     END IF

  ELSE rnd_store
     ! Re-use stored rndmat
     if (mype == 0 .AND. screen > 0) &
          write (*,'(a, 5x, a)') 'PDAF', '--- Use stored random rotation matrix'
     rndmat = rndmat_save

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  stored rndmat', rndmat
  END IF rnd_store

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

  CALL PDAF_timeit(4, 'old')


! *************************************
! *** Perform analysis - first loop ***
! *************************************

  CALL PDAF_timeit(6, 'new')

  ! Initialize counters for statistics on local observations
  CALL PDAF_init_local_obsstats()

!$OMP PARALLEL default(shared) private(dim_l, dim_obs_l, obs_l, HX_l, HXbar_l, resid_l, ens_l) &
!$OMP private(state_l, stateinc_l, Uinv_l, forget_ana_l)

  forget_ana_l = forget_ana

  ! Allocate ensemble transform matrix
  ALLOCATE(Uinv_l(dim_ens, dim_ens))
  IF (allocflag == 0) CALL PDAF_memcount(3, 'r', dim_ens**2)
  Uinv_l = 0.0

  ! initialize number of small singular values
  cnt_small_svals = 0

  IF (debug>0 .and. n_domains_p>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- Enter first local analysis loop'

!$OMP BARRIER
!$OMP DO firstprivate(cnt_maxlag) lastprivate(cnt_maxlag) schedule(runtime)
  localanalysis: DO domain_p = 1, n_domains_p

     ! Set flag that we are in the local analysis loop
     inloop = .true.

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- First local analysis for domain_p', domain_p
        WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_dim_l'
     END IF

     ! local state dimension
     CALL PDAF_timeit(45, 'new')
     CALL U_init_dim_l(step, domain_p, dim_l)
     CALL PDAF_timeit(45, 'old')

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  dim_l', dim_l
        WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_dim_obs_l'
     END IF

     ! Get observation dimension for local domain
     CALL PDAF_timeit(9, 'new')
     dim_obs_l = 0
     CALL U_init_dim_obs_l(domain_p, step, dim_obs_f, dim_obs_l)
     CALL PDAF_timeit(9, 'old')

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  dim_obs_l', dim_obs_l

     CALL PDAF_timeit(51, 'new')

     ! Gather statistical information on local observations
     CALL PDAF_incr_local_obsstats(dim_obs_l)
     
     ! Allocate arrays for local analysis domain
     ALLOCATE(ens_l(dim_l, dim_ens))
     ALLOCATE(state_l(dim_l))
     ALLOCATE(stateinc_l(dim_l))
     CALL PDAF_timeit(51, 'old')

     CALL PDAF_timeit(15, 'new')

     ! state ensemble and mean state on current analysis domain
     DO member = 1, dim_ens
        ! Store member index to make it accessible with PDAF_get_memberid
        member_save = member

        IF (debug>0) then
           WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update -- call g2l_state for ensemble member', member
           if (member==1) &
                WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update --    Note: if ens_l is incorrect check user-defined indices in g2l_state!'
        END IF

        CALL U_g2l_state(step, domain_p, dim_p, ens_p(:, member), dim_l, &
             ens_l(:, member))

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  ens_l', ens_l(:,member)
     END DO

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, &
          'PDAF_lknetf_update -- call g2l_state for ensemble mean'

     member_save = 0
     CALL U_g2l_state(step, domain_p, dim_p, state_p, dim_l, &
          state_l)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  meanens_l', state_l

     CALL PDAF_timeit(15, 'old')

     IF (dim_obs_l > 0) THEN

        ALLOCATE(HX_l(dim_obs_l, dim_ens))
        ALLOCATE(HXbar_l(dim_obs_l))
        ALLOCATE(obs_l(dim_obs_l))

        ! Restrict full observation to local observation
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- call g2l_obs', dim_ens, 'times'

        CALL PDAF_timeit(46, 'new')
        DO member = 1, dim_ens
           obs_member = member
           CALL U_g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, HX_f(:,member), HX_l(:,member))
        END DO

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call g2l_obs for mean'

        obs_member = 0
        CALL U_g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, HXbar_f, HXbar_l)
        CALL PDAF_timeit(46, 'old')

        ! get local observation vector
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_obs_l'

        CALL PDAF_timeit(47, 'new')
        CALL U_init_obs_l(domain_p, step, dim_obs_l, obs_l)
        CALL PDAF_timeit(47, 'old')

        ! Omit observations with too high innovation
        IF (omi_omit_obs)  THEN
           CALL PDAF_timeit(51, 'new')

           ALLOCATE(resid_l(dim_obs_l))

           resid_l = obs_l - HXbar_l
           CALL PDAFomi_omit_by_inno_l_cb(domain_p, dim_obs_l, resid_l, obs_l)

           DEALLOCATE(resid_l)

           CALL PDAF_timeit(51, 'old')
        END IF

        CALL PDAF_timeit(7, 'new')


        ! *** 2-step LKNETF analysis - STEP 1 ***

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

        CALL PDAF_timeit(53,'new')
        CALL PDAF_lknetf_compute_gamma(domain_p, step, dim_obs_l, dim_ens, &
             HX_l, HXbar_l, obs_l, type_hyb, hyb_g, hyb_k, &
             gamma(domain_p), n_eff(domain_p), skewness(domain_p), kurtosis(domain_p), &
             U_likelihood_l, screen, flag)
        CALL PDAF_timeit(53,'old')

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

        IF (subtype == 0) THEN

           ! 2-step LKNETF with NETF before LETKF 
           CALL PDAF_lknetf_ana_lnetf(domain_p, step, dim_l, dim_obs_l, &
                dim_ens, ens_l, HX_l, rndmat, obs_l, U_likelihood_hyb_l, &
                cnt_small_svals, n_eff_b(domain_p), gamma(domain_p), screen, flag)
        ELSE IF (subtype == 1) THEN

           ! 2-step LKNETF with LETKF before NETF
           CALL PDAF_lknetf_ana_letkfT(domain_p, step, dim_l, dim_obs_l, &
                dim_ens, state_l, Uinv_l, ens_l, HX_l, &
                HXbar_l, stateinc_l, rndmat, forget_ana_l, &
                obs_l, U_prodRinvA_hyb_l, U_init_obsvar_l, U_init_n_domains_p, &
                gamma(domain_p), screen, incremental, type_forget, flag)
        END IF

        DEALLOCATE(obs_l, HX_l, HXbar_l)

        CALL PDAF_timeit(7, 'old')

     ELSE

        ! UNOBSERVED DOMAIN

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- dim_obs_l = 0; omit call to local analysis function'

        IF (type_forget==1) THEN
           ! prior inflation NOT on unobserved domains - take it back!
           invforget = 1.0/forget
           CALL PDAF_inflate_ens(dim_l, dim_ens, state_p, ens_l, invforget)
        ENDIF

     END IF

     CALL PDAF_timeit(16, 'new')
 
     ! re-initialize full state ensemble on PE and mean state from local domain
     DO member = 1, dim_ens
        IF (debug>0) then
           WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update -- call l2g_state for ensemble member', member
           WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  ens_l', ens_l(:,member)
        END IF

        member_save = member
        CALL U_l2g_state(step, domain_p, dim_l, ens_l(:, member), dim_p, ens_p(:,member))
     END DO
    
     ! Initialize global state increment
!      IF (incremental == 1) THEN
!         IF (debug>0) THEN
!            WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- init gobal state increment'
!            WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  stateinc_l', stateinc_l
!         END IF
!
!         member_save = member
!         CALL U_l2g_state(step, domain_p, dim_l, stateinc_l, dim_p, state_inc_p)
!      END IF

     CALL PDAF_timeit(16, 'old')

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

     ! *** Perform smoothing of past ensembles ***
!     CALL PDAF_smoother_local(domain_p, step, dim_p, dim_l, dim_ens, &
!          dim_lag, Uinv_l, ens_l, sens_p, cnt_maxlag, &
!          U_g2l_state, U_l2g_state, forget_ana, screen)

     CALL PDAF_timeit(17, 'old')

     ! clean up
     DEALLOCATE(ens_l, state_l, stateinc_l)
     CALL PDAF_timeit(51, 'old')

  END DO localanalysis

  IF (debug>0 .and. n_domains_p>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- End of first local analysis loop'

  ! Set flag that we are not in the local analysis loop
  inloop = .false.


! ***************************************
! *** Prepare second step of analysis ***
! ***************************************

  CALL PDAF_timeit(4, 'new')

  ! HX = [Hx_1 ... Hx_(r+1)] for full DIM_OBS_F region on PE-local domain
  CALL PDAF_timeit(44, 'new')

  IF (debug>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call obs_op', dim_ens, 'times'

  ENSA: DO member = 1,dim_ens
     ! Store member index to make it accessible with PDAF_get_obsmemberid
     obs_member = member

     ! Call observation operator
     CALL U_obs_op(step, dim_p, dim_obs_f, ens_p(:, member), HX_f(:, member))
  END DO ENSA
  CALL PDAF_timeit(44, 'old')

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

  ! *** Compute mean forecast state
  state_p = 0.0
  invdimens = 1.0 / REAL(dim_ens)
  DO member = 1, dim_ens
     DO row = 1, dim_p
        state_p(row) = state_p(row) + invdimens * ens_p(row, member)
     END DO
  END DO

  CALL PDAF_timeit(11, 'old')

  ! *** Compute mean state of ensemble on PE-local observation space 
  HXbar_f = 0.0
  invdimens = 1.0 / REAL(dim_ens)
  DO member = 1, dim_ens
     DO row = 1, dim_obs_f
        HXbar_f(row) = HXbar_f(row) + invdimens * HX_f(row, member)
     END DO
  END DO
  CALL PDAF_timeit(51, 'old')

  CALL PDAF_timeit(4, 'old')


! **************************************
! *** Perform analysis - second loop ***
! **************************************

  CALL PDAF_timeit(6, 'new')

  forget_ana_l = forget_ana

  ! Allocate ensemble transform matrix
  Uinv_l = 0.0

  ! initialize number of small singular values
  cnt_small_svals = 0

  IF (debug>0 .and. n_domains_p>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- Enter second local analysis loop'

!$OMP BARRIER
!$OMP DO firstprivate(cnt_maxlag) lastprivate(cnt_maxlag) schedule(runtime)
  localanalysisA: DO domain_p = 1, n_domains_p

     ! Set flag that we are in the local analysis loop
     inloop = .true.

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- Second local analysis for domain_p', domain_p
        WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_dim_l'
     END IF

     ! local state dimension
     CALL PDAF_timeit(45, 'new')
     CALL U_init_dim_l(step, domain_p, dim_l)
     CALL PDAF_timeit(45, 'old')

     IF (debug>0) THEN
        WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  dim_l', dim_l
        WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_dim_obs_l'
     END IF

     ! Get observation dimension for local domain
     CALL PDAF_timeit(9, 'new')
     dim_obs_l = 0
     CALL U_init_dim_obs_l(domain_p, step, dim_obs_f, dim_obs_l)
     CALL PDAF_timeit(9, 'old')

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  dim_obs_l', dim_obs_l

     CALL PDAF_timeit(51, 'new')
     ! Allocate arrays for local analysis domain
     ALLOCATE(ens_l(dim_l, dim_ens))
     ALLOCATE(state_l(dim_l))
     ALLOCATE(stateinc_l(dim_l))
     CALL PDAF_timeit(51, 'old')

     CALL PDAF_timeit(15, 'new')

     ! state ensemble and mean state on current analysis domain
     DO member = 1, dim_ens
        IF (debug>0) then
           WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update -- call g2l_state for ensemble member', member
           if (member==1) &
                WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update --    Note: if ens_l is incorrect check user-defined indices in g2l_state!'
        END IF

        member_save = member
        CALL U_g2l_state(step, domain_p, dim_p, ens_p(:, member), dim_l, &
             ens_l(:, member))

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  ens_l', ens_l(:,member)
     END DO

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug: ', debug, &
          'PDAF_lknetf_update -- call g2l_state for ensemble mean'

     member_save = 0
     CALL U_g2l_state(step, domain_p, dim_p, state_p, dim_l, &
          state_l)

     IF (debug>0) &
          WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  meanens_l', state_l

     CALL PDAF_timeit(15, 'old')

     IF (dim_obs_l > 0) THEN

        ! get local observation vector
        ALLOCATE(HX_l(dim_obs_l, dim_ens))
        ALLOCATE(HXbar_l(dim_obs_l))
        ALLOCATE(obs_l(dim_obs_l))

        ! Restrict full observation to local observation
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- call g2l_obs', dim_ens, 'times'

        CALL PDAF_timeit(46, 'new')
        DO member = 1, dim_ens
           obs_member = member
           CALL U_g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, HX_f(:,member), HX_l(:,member))
        END DO

        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call g2l_obs for mean'

        obs_member = 0
        CALL U_g2l_obs(domain_p, step, dim_obs_f, dim_obs_l, HXbar_f, HXbar_l)
        CALL PDAF_timeit(46, 'old')

        ! get local observation vector
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- call init_obs_l'

        CALL PDAF_timeit(47, 'new')
        CALL U_init_obs_l(domain_p, step, dim_obs_l, obs_l)
        CALL PDAF_timeit(47, 'new')

        CALL PDAF_timeit(7, 'new')


        ! *** 2-step LKNETF analysis - STEP 2 ***

        IF (subtype == 0) THEN

           ! 2-step LKNETF with NETF before LETKF 
           CALL PDAF_lknetf_ana_letkfT(domain_p, step, dim_l, dim_obs_l, &
                dim_ens, state_l, Uinv_l, ens_l, HX_l, &
                HXbar_l, stateinc_l, rndmat, forget_ana_l, &
                obs_l, U_prodRinvA_hyb_l, U_init_obsvar_l, U_init_n_domains_p, &
                gamma(domain_p), screen, incremental, type_forget, flag)
        ELSE IF (subtype == 1) THEN
           CALL PDAF_lknetf_ana_lnetf(domain_p, step, dim_l, dim_obs_l, &
                dim_ens, ens_l, HX_l, rndmat, obs_l, U_likelihood_hyb_l, &
                cnt_small_svals, n_eff_b(domain_p), gamma(domain_p), screen, flag)
        END IF

        IF (type_forget==3) THEN
           IF (debug>0) &
                WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- apply inflation to analysis ensemble'

           ! Apply forgetting factor to posterior ensemble - only observed domains
           CALL PDAF_timeit(14, 'new')
           CALL PDAF_inflate_ens(dim_l, dim_ens, state_l, ens_l, forget)
           CALL PDAF_timeit(14, 'old')
        ENDIF

        CALL PDAF_timeit(7, 'old')

        DEALLOCATE(obs_l, HX_l, HXbar_l)

     ELSE
        IF (debug>0) &
             WRITE (*,*) '++ PDAF-debug: ', debug, &
             'PDAF_lknetf_update -- dim_obs_l = 0; omit call to local analysis function'
     END IF

     CALL PDAF_timeit(16, 'new')
 
     ! re-initialize full state ensemble on PE and mean state from local domain
     DO member = 1, dim_ens
        IF (debug>0) then
           WRITE (*,*) '++ PDAF-debug: ', debug, &
                'PDAF_lknetf_update -- call l2g_state for ensemble member', member
           WRITE (*,*) '++ PDAF-debug PDAF_lknetf_update:', debug, '  ens_l', ens_l(:,member)
        END IF

        member_save = member
        CALL U_l2g_state(step, domain_p, dim_l, ens_l(:, member), dim_p, ens_p(:,member))
     END DO
    
     ! Initialize global state increment
!      IF (incremental == 1) THEN
!         member_save = member
!         CALL U_l2g_state(step, domain_p, dim_l, stateinc_l, dim_p, state_inc_p)
!      END IF

     CALL PDAF_timeit(16, 'old')

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

     ! *** Perform smoothing of past ensembles ***
!     CALL PDAF_smoother_local(domain_p, step, dim_p, dim_l, dim_ens, &
!          dim_lag, Uinv_l, ens_l, sens_p, cnt_maxlag, &
!          U_g2l_state, U_l2g_state, forget_ana, screen)

     CALL PDAF_timeit(17, 'old')

     ! clean up
     DEALLOCATE(ens_l, state_l, stateinc_l)
     CALL PDAF_timeit(51, 'old')

  END DO localanalysisA

  IF (debug>0 .and. n_domains_p>0) &
       WRITE (*,*) '++ PDAF-debug: ', debug, 'PDAF_lknetf_update -- End of second local analysis loop'

  ! Set flag that we are not in the local analysis loop
  inloop = .false.

  CALL PDAF_timeit(51, 'new')

  IF (type_forget==2) THEN
     ! Apply forgetting factor to posterior ensemble - all domains
     CALL PDAF_timeit(14, 'new')
     CALL PDAF_inflate_ens(dim_p, dim_ens, state_p, ens_p, forget)
     CALL PDAF_timeit(14, 'old')
  ENDIF

!$OMP CRITICAL
  ! Set Uinv - required for subtype=3
  Uinv = Uinv_l
!$OMP END CRITICAL

  DEALLOCATE(Uinv_l)
!$OMP END PARALLEL

  CALL PDAF_timeit(6, 'old')

  CALL PDAF_timeit(3, 'old')

  ! Initialize mask array for computing global effective ensemble size
  MASK = (n_eff_b > 0)

  ! *** Print statistics for local analysis to the screen ***
  CALL PDAF_print_local_obsstats(screen, n_domains_with_obs_p)

  IF (npes_filter>1) THEN
     ! Number globally observed domains
     CALL MPI_Reduce(n_domains_with_obs_p, n_domains_with_obs, 1, MPI_INTEGER, MPI_SUM, &
          0, COMM_filter, MPIerr)

     ! Min/max effective sample sizes
     max_n_eff_l = MAXVAL(n_eff_b)
     CALL MPI_Reduce(max_n_eff_l, max_n_eff, 1, MPI_REALTYPE, MPI_MAX, &
          0, COMM_filter, MPIerr)
     min_n_eff_l = MINVAL(n_eff_b, MASK)
     CALL MPI_Reduce(min_n_eff_l, min_n_eff, 1, MPI_REALTYPE, MPI_MIN, &
          0, COMM_filter, MPIerr)
     sum_n_eff_l = SUM(n_eff_b)
     CALL MPI_Reduce(sum_n_eff_l, mean_n_eff, 1, MPI_REALTYPE, MPI_SUM, &
          0, COMM_filter, MPIerr)
     mean_n_eff = mean_n_eff / REAL(n_domains_with_obs)

     ! Min/max hybrid weight
     max_gamma_l = MAXVAL(gamma)
     CALL MPI_Reduce(max_gamma_l, max_gamma, 1, MPI_REALTYPE, MPI_MAX, &
          0, COMM_filter, MPIerr)
     min_gamma_l = MINVAL(gamma, MASK)
     CALL MPI_Reduce(min_gamma_l, min_gamma, 1, MPI_REALTYPE, MPI_MIN, &
          0, COMM_filter, MPIerr)
     sum_gamma_l = SUM(gamma)
     CALL MPI_Reduce(sum_gamma_l, mean_gamma, 1, MPI_REALTYPE, MPI_SUM, &
          0, COMM_filter, MPIerr)
     mean_gamma = mean_gamma / REAL(n_domains_with_obs)

     ! Min/max skewness and kurtosis
     max_stats_l(1) = MAXVAL(skewness)/SQRT(REAL(dim_ens))
     max_stats_l(2) = MAXVAL(kurtosis)/REAL(dim_ens)
     CALL MPI_Reduce(max_stats_l, max_stats, 2, MPI_REALTYPE, MPI_MAX, &
          0, COMM_filter, MPIerr)

     min_stats_l(1) = MINVAL(skewness, MASK)/SQRT(REAL(dim_ens))
     min_stats_l(2) = MINVAL(kurtosis, MASK)/REAL(dim_ens)
     CALL MPI_Reduce(min_stats_l, min_stats, 2, MPI_REALTYPE, MPI_MIN, &
          0, COMM_filter, MPIerr)

     sum_stats_l(1) = SUM(skewness)/SQRT(REAL(dim_ens))
     sum_stats_l(2) = SUM(kurtosis)/REAL(dim_ens)
     CALL MPI_Reduce(sum_stats_l, mean_stats, 2, MPI_REALTYPE, MPI_SUM, &
          0, COMM_filter, MPIerr)
     mean_stats = mean_stats / REAL(n_domains_with_obs)

  ELSE
     ! Min/max effective ensemble sizes
     max_n_eff = MAXVAL(n_eff_b)
     min_n_eff = MINVAL(n_eff_b, MASK)
     mean_n_eff = SUM(n_eff_b) / n_domains_with_obs_p

     ! Min/max effective ensemble sizes
     max_gamma = MAXVAL(gamma)
     min_gamma = MINVAL(gamma, MASK)
     mean_gamma = SUM(gamma) / n_domains_with_obs_p

     ! Min/max skewness and kurtosis
     max_stats(1) = MAXVAL(skewness)/SQRT(REAL(dim_ens))
     max_stats(2) = MAXVAL(kurtosis)/REAL(dim_ens)
     min_stats(1) = MINVAL(skewness, MASK)/SQRT(REAL(dim_ens))
     min_stats(2) = MINVAL(kurtosis, MASK)/REAL(dim_ens)
     mean_stats(1) = SUM(skewness)/SQRT(REAL(dim_ens)) / REAL(n_domains_with_obs_p)
     mean_stats(2) = SUM(kurtosis)/REAL(dim_ens)/ REAL(n_domains_with_obs_p)
  END IF

  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 8x, a, 23x, 3f10.3)') &
         'PDAF', 'Minimal/Maximal/Mean N_eff:', &
         min_n_eff, max_n_eff, mean_n_eff
     WRITE (*, '(a, 8x, a, 15x, 3f10.3)') &
         'PDAF', 'Minimal/Maximal/Mean hybrid weight:', &
         min_gamma, max_gamma, mean_gamma
     WRITE (*, '(a, 8x, a, 1x, 3f10.3)') &
         'PDAF', 'Minimal/Maximal/Mean abs. skewness/SQRT(dim_ens):', &
         min_stats(1), max_stats(1), mean_stats(1)
     WRITE (*, '(a, 8x, a, 7x, 3f10.3)') &
         'PDAF', 'Minimal/Maximal/Mean abs. kurtosis/dim_ens:', &
         min_stats(2), max_stats(2), mean_stats(2)

     IF (screen > 1) THEN
        WRITE (*, '(a, 5x, a, F10.3, 1x, a)') &
             'PDAF', '--- analysis/re-init duration:', PDAF_time_temp(3), 's'
     END IF
  END IF
  CALL PDAF_timeit(51, 'old')

! *** Clean up from local analysis update ***
  DEALLOCATE(HX_f, HXbar_f, rndmat)
#else
  WRITE (*,'(/5x,a/)') &
       '!!! PDAF WARNING: ANALYSIS STEP IS DEACTIVATED BY PDAF_NO_UPDATE !!!'
#endif

! *** Poststep for analysis ensemble ***
  CALL PDAF_timeit(5, 'new')
  IF (mype == 0 .AND. screen > 0) THEN
     WRITE (*, '(a, 5x, a)') 'PDAF', 'Call pre-post routine after analysis step'
  ENDIF
  CALL U_prepoststep(step, dim_p, dim_ens, dim_ens_l, dim_obs_f, &
       state_p, Uinv, ens_p, flag)
  CALL PDAF_timeit(5, 'old')
  
  IF (mype == 0 .AND. screen > 0) THEN
     IF (screen > 1) THEN
        WRITE (*, '(a, 5x, a, F10.3, 1x, a)') &
             'PDAF', '--- duration of poststep:', PDAF_time_temp(5), 's'
     END IF
     WRITE (*, '(a, 55a)') 'PDAF Forecast ', ('-', i = 1, 55)
  END IF


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

  DEALLOCATE(n_eff, n_eff_b, MASK, gamma)
  DEALLOCATE(skewness, kurtosis)

  IF (allocflag == 0) allocflag = 1

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

END SUBROUTINE PDAF_lknetf_step_update