Merge pull request #794 from NCAR/qceff_inflation_fix

Qceff inflation fix; Explicitly handling BNRHF transform failures.

Author: hkershaw-brown

CHANGELOG.rst

@@ -22,6 +22,11 @@ individual files.
 
 The changes are now listed with the most recent at the top.
 
+**January 30 2025 :: Bug-fix: Explicitly handle BNRHF transform failures. Tag v11.10.1**
+
+- Probit transform failure is caught and an error code is returned 
+- filter_mod and assim_tools_mod skip variables that fail the transform
+
 **January 23 2025 :: DART_LAB QCEFF. Tag v11.10.0**
 
 - Updated DART_LAB to include QCEFF

assimilation_code/modules/assimilation/assim_tools_mod.f90

@@ -346,7 +346,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
 integer(i8), allocatable :: my_state_indx(:)
 integer(i8), allocatable :: my_obs_indx(:)
 
-integer :: my_num_obs, i, j, owner, owners_index, my_num_state
+integer :: my_num_obs, i, j, owner, owners_index, my_num_state, ierr
 integer :: obs_mean_index, obs_var_index
 integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group
 integer :: num_close_obs, obs_index, num_close_states
@@ -373,6 +373,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
 type(obs_def_type)   :: obs_def
 type(time_type)      :: obs_time
 
+logical, allocatable :: obs_probit_trans_ok(:), state_probit_trans_ok(:)
 logical :: allow_missing_in_state
 logical :: local_single_ss_inflate
 logical :: local_varying_ss_inflate
@@ -395,13 +396,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
          my_obs_indx(        obs_ens_handle%my_num_vars), &
          my_obs_kind(        obs_ens_handle%my_num_vars), &
          my_obs_type(        obs_ens_handle%my_num_vars), &
-         my_obs_loc(         obs_ens_handle%my_num_vars))
+         my_obs_loc(         obs_ens_handle%my_num_vars), &
+         obs_probit_trans_ok(obs_ens_handle%my_num_vars))
 
 allocate(close_state_dist(     ens_handle%my_num_vars), &
          close_state_ind(      ens_handle%my_num_vars), &
          my_state_indx(        ens_handle%my_num_vars), &
          my_state_kind(        ens_handle%my_num_vars), &
-         my_state_loc(         ens_handle%my_num_vars))
+         my_state_loc(         ens_handle%my_num_vars), &
+         state_probit_trans_ok(ens_handle%my_num_vars))
 ! end alloc
 
 ! Initialize assim_tools_module if needed
@@ -503,8 +506,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
    ! Convert all my state variables to appropriate probit space
    call transform_to_probit(ens_size, ens_handle%copies(1:ens_size, i), dist_for_state, &
       state_dist_params(i), probit_ens, .false., &
-      bounded_below, bounded_above, lower_bound, upper_bound)
-   ens_handle%copies(1:ens_size, i) = probit_ens
+      bounded_below, bounded_above, lower_bound, upper_bound, ierr)
+   state_probit_trans_ok(i) = (ierr == 0)
+   if(state_probit_trans_ok(i)) ens_handle%copies(1:ens_size, i) = probit_ens
 end do
 
 !> optionally convert all state location verticals
@@ -541,8 +545,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
       ! Convert all my obs (extended state) variables to appropriate probit space
       call transform_to_probit(ens_size, obs_ens_handle%copies(1:ens_size, i), dist_for_obs, &
          obs_dist_params(i), probit_ens, .false., &
-         bounded_below, bounded_above, lower_bound, upper_bound)
-      obs_ens_handle%copies(1:ens_size, i) = probit_ens
+         bounded_below, bounded_above, lower_bound, upper_bound, ierr)
+      obs_probit_trans_ok(i) = (ierr == 0)
+      if(obs_probit_trans_ok(i)) obs_ens_handle%copies(1:ens_size, i) = probit_ens
    endif
 end do
 
@@ -635,9 +640,11 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
          orig_obs_prior_var  = obs_ens_handle%copies(OBS_PRIOR_VAR_START:  &
             OBS_PRIOR_VAR_END, owners_index)
 
-         ! If QC is okay, convert this observation ensemble from probit to regular space
-         call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , &
-            obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index))
+         ! Convert this observation ensemble from probit back to regular space if to_probit was successful
+         if(obs_probit_trans_ok(owners_index)) then 
+            call transform_from_probit(ens_size, obs_ens_handle%copies(1:ens_size, owners_index) , &
+               obs_dist_params(owners_index), obs_ens_handle%copies(1:ens_size, owners_index))
+         endif
 
          obs_prior = obs_ens_handle%copies(1:ens_size, owners_index)
       endif IF_QC_IS_OKAY
@@ -698,10 +705,14 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
       ! Convert the prior and posterior for this observation to probit space
       call transform_to_probit(grp_size, obs_prior(grp_bot:grp_top), dist_for_obs, &
          temp_dist_params, probit_obs_prior(grp_bot:grp_top), .false., &
-         bounded_below, bounded_above, lower_bound, upper_bound)
+         bounded_below, bounded_above, lower_bound, upper_bound, ierr)
+      ! If probit transform fails for any group, this observation cannot be assimimlated
+      if(ierr /= 0) cycle SEQUENTIAL_OBS
       call transform_to_probit(grp_size, obs_post(grp_bot:grp_top), dist_for_obs, &
          temp_dist_params, probit_obs_post(grp_bot:grp_top), .true., &
-         bounded_below, bounded_above, lower_bound, upper_bound)
+         bounded_below, bounded_above, lower_bound, upper_bound, ierr)
+      ! If probit transform fails for any group, this observation cannot be assimimlated
+      if(ierr /= 0) cycle SEQUENTIAL_OBS
       ! Free up the storage used for this transform
       call deallocate_distribution_params(temp_dist_params)
 
@@ -761,6 +772,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
    STATE_UPDATE: do j = 1, num_close_states
       state_index = close_state_ind(j)
 
+      ! If transform to probit failed for this state variable, don't update it
+      if(.not. state_probit_trans_ok(state_index)) cycle STATE_UPDATE
+
       if ( allow_missing_in_state ) then
          ! Don't allow update of state ensemble with any missing values
          if (any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)) cycle STATE_UPDATE
@@ -796,6 +810,9 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
       OBS_UPDATE: do j = 1, num_close_obs
          obs_index = close_obs_ind(j)
 
+         ! If transform to probit failed for this observed variable, don't update it
+         if(.not. obs_probit_trans_ok(obs_index)) cycle OBS_UPDATE
+
          ! Only have to update obs that have not yet been used
          if(my_obs_indx(obs_index) > i) then
 
@@ -820,7 +837,7 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
 
 ! Do the inverse probit transform for state variables
 call transform_all_from_probit(ens_size, ens_handle%my_num_vars, ens_handle%copies, &
-   state_dist_params, ens_handle%copies)
+   state_dist_params, ens_handle%copies, state_probit_trans_ok)
 
 ! Every pe needs to get the current my_inflate and my_inflate_sd back
 if(local_single_ss_inflate) then
@@ -876,13 +893,15 @@ subroutine filter_assim(ens_handle, obs_ens_handle, obs_seq, keys,           &
            my_obs_type,         &
            close_obs_ind,       &
            vstatus,             &
-           my_obs_loc)
+           my_obs_loc,          &
+           state_probit_trans_ok)
 
 deallocate(close_state_dist,      &
            my_state_indx,         &
            close_state_ind,       &
            my_state_kind,         &
-           my_state_loc)
+           my_state_loc,          &
+           obs_probit_trans_ok)
 
 ! end dealloc
 

assimilation_code/modules/assimilation/filter_mod.f90

@@ -1567,7 +1567,7 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C
 real(r8) :: probit_ens(ens_size), probit_ens_mean
 logical  :: bounded_below, bounded_above
 real(r8) :: lower_bound,   upper_bound
-integer  :: dist_type
+integer  :: dist_type, ierr
 
 ! Inflate each group separately;  Divide ensemble into num_groups groups
 grp_size = ens_size / num_groups
@@ -1609,16 +1609,20 @@ subroutine filter_ensemble_inflate(ens_handle, inflate_copy, inflate, ENS_MEAN_C
 
          call transform_to_probit(grp_size, ens_handle%copies(grp_bot:grp_top, j), &
             dist_type, dist_params, probit_ens(1:grp_size), .false., &
-               bounded_below, bounded_above, lower_bound, upper_bound)
-
-         ! Compute the ensemble mean in transformed space
-         probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size
-         ! Inflate in probit space
-         call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, &
-            ens_handle%copies(inflate_copy, j))
-         ! Transform back from probit space
-         call transform_from_probit(grp_size, probit_ens(1:grp_size), &
-            dist_params, ens_handle%copies(grp_bot:grp_top, j))
+               bounded_below, bounded_above, lower_bound, upper_bound, ierr)
+
+         ! Only inflate if successful return from the probit transform
+         if(ierr == 0) then
+            ! Compute the ensemble mean in transformed space
+            probit_ens_mean = sum(probit_ens(1:grp_size)) / grp_size
+            ! Inflate in probit space
+            call inflate_ens(inflate, probit_ens(1:grp_size), probit_ens_mean, &
+               ens_handle%copies(inflate_copy, j))
+            ! Transform back from probit space
+            call transform_from_probit(grp_size, probit_ens(1:grp_size), &
+               dist_params, ens_handle%copies(grp_bot:grp_top, j))
+         endif
+
       end do
    endif
 end do