[Dart-dev] [4959] DART/trunk/filter: Changes in how we set the mean and stddev when there are failed

nancy at ucar.edu nancy at ucar.edu
Wed Jun 8 08:55:02 MDT 2011


Revision: 4959
Author:   nancy
Date:     2011-06-08 08:55:02 -0600 (Wed, 08 Jun 2011)
Log Message:
-----------
Changes in how we set the mean and stddev when there are failed
forward operators.  The older code set the values right before
they were written out.  This code sets the values inside the obs
loop right when they're computed, and in addition has a section
to set the posterior values if any of the posterior forward operators
failed, even if the prior qc is bad (the particular case where this
is possible is when the prior fails the outlier threshold test, but
all the prior forward operators worked).   Also make the same changes
to the radar folded version of filter to keep it in sync.

Change the names of OBS_PRIOR_MEAN_... and OBS_PRIOR_VAR_...
to just OBS_MEAN_... and OBS_VAR_... because they are used both
when computing the prior and posterior values.  these are internal
to the filter_main() subroutine and so are not visible outside the
subroutines in which they are used.  just for a bit more clarity.
(there's already a separate flag passed down to the obs diag routine 
to say whether we're computing the prior or posterior.)

Modified Paths:
--------------
    DART/trunk/filter/filter.dopplerfold.f90
    DART/trunk/filter/filter.f90

-------------- next part --------------
Modified: DART/trunk/filter/filter.dopplerfold.f90
===================================================================
--- DART/trunk/filter/filter.dopplerfold.f90	2011-06-07 22:54:21 UTC (rev 4958)
+++ DART/trunk/filter/filter.dopplerfold.f90	2011-06-08 14:55:02 UTC (rev 4959)
@@ -181,8 +181,8 @@
 integer                 :: ENS_MEAN_COPY, ENS_SD_COPY, PRIOR_INF_COPY, PRIOR_INF_SD_COPY
 integer                 :: POST_INF_COPY, POST_INF_SD_COPY
 integer                 :: OBS_VAL_COPY, OBS_ERR_VAR_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY
-integer                 :: OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END
-integer                 :: OBS_PRIOR_VAR_START, OBS_PRIOR_VAR_END, TOTAL_OBS_COPIES
+integer                 :: OBS_MEAN_START, OBS_MEAN_END
+integer                 :: OBS_VAR_START, OBS_VAR_END, TOTAL_OBS_COPIES
 integer                 :: input_qc_index, DART_qc_index
 integer                 :: mean_owner, mean_owners_index
 
@@ -245,10 +245,10 @@
 OBS_VAL_COPY         = ens_size + 2
 OBS_KEY_COPY         = ens_size + 3
 OBS_GLOBAL_QC_COPY   = ens_size + 4
-OBS_PRIOR_MEAN_START = ens_size + 5
-OBS_PRIOR_MEAN_END   = OBS_PRIOR_MEAN_START + num_groups - 1
-OBS_PRIOR_VAR_START  = OBS_PRIOR_MEAN_START + num_groups
-OBS_PRIOR_VAR_END    = OBS_PRIOR_VAR_START + num_groups - 1
+OBS_MEAN_START       = ens_size + 5
+OBS_MEAN_END         = OBS_MEAN_START + num_groups - 1
+OBS_VAR_START        = OBS_MEAN_START + num_groups
+OBS_VAR_END          = OBS_VAR_START + num_groups - 1
 
 ! Can't output more ensemble members than exist
 if(num_output_state_members > ens_size) num_output_state_members = ens_size
@@ -547,7 +547,7 @@
       seq, keys, PRIOR_DIAG, num_output_obs_members, in_obs_copy+1, &
       obs_val_index, OBS_KEY_COPY, &                                 ! new
       prior_obs_mean_index, prior_obs_spread_index, num_obs_in_set, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, &
+      OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
    call trace_message('After  observation space diagnostics')
   
@@ -564,8 +564,8 @@
       ens_size, num_groups, obs_val_index, prior_inflate, &
       ENS_MEAN_COPY, ENS_SD_COPY, &
       PRIOR_INF_COPY, PRIOR_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-      OBS_PRIOR_VAR_END, inflate_only = .false.)
+      OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+      OBS_VAR_END, inflate_only = .false.)
 
    call timestamp_message('After  observation assimilation')
    call     trace_message('After  observation assimilation')
@@ -582,8 +582,8 @@
       call smoother_assim(obs_ens_handle, seq, keys, ens_size, num_groups, &
          obs_val_index, ENS_MEAN_COPY, ENS_SD_COPY, &
          PRIOR_INF_COPY, PRIOR_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-         OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-         OBS_PRIOR_VAR_END)
+         OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+         OBS_VAR_END)
       call timestamp_message('After  smoother assimilation')
       call     trace_message('After  smoother assimilation')
    endif
@@ -656,7 +656,7 @@
       seq, keys, POSTERIOR_DIAG, num_output_obs_members, in_obs_copy+2, &
       obs_val_index, OBS_KEY_COPY, &                             ! new
       posterior_obs_mean_index, posterior_obs_spread_index, num_obs_in_set, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, &
+      OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
    
    call trace_message('After  posterior obs space diagnostics')
@@ -692,8 +692,8 @@
          call filter_assim(ens_handle, obs_ens_handle, seq, keys, ens_size, num_groups, &
             obs_val_index, post_inflate, ENS_MEAN_COPY, ENS_SD_COPY, &
             POST_INF_COPY, POST_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-            OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-            OBS_PRIOR_VAR_END, inflate_only = .true.)
+            OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+            OBS_VAR_END, inflate_only = .true.)
 
          call all_copies_to_all_vars(ens_handle)
 
@@ -1315,7 +1315,7 @@
    seq, keys, prior_post, num_output_members, members_index, &
    obs_val_index, OBS_KEY_COPY, &
    ens_mean_index, ens_spread_index, num_obs_in_set, &
-   OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
+   OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
    OBS_ERR_VAR_COPY, DART_qc_index)
 
 ! Do prior observation space diagnostics on the set of obs corresponding to keys
@@ -1329,17 +1329,17 @@
 integer,                 intent(in)    :: OBS_KEY_COPY
 integer,                 intent(in)    :: ens_mean_index, ens_spread_index
 type(obs_sequence_type), intent(inout) :: seq
-integer,                 intent(in)    :: OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START
+integer,                 intent(in)    :: OBS_MEAN_START, OBS_VAR_START
 integer,                 intent(in)    :: OBS_GLOBAL_QC_COPY, OBS_VAL_COPY
 integer,                 intent(in)    :: OBS_ERR_VAR_COPY, DART_qc_index
 
 integer               :: j, k, ens_offset, forward_min, forward_max
 integer               :: forward_unit, ivalue, low_qc_limit, high_qc_limit
 real(r8)              :: error, diff_sd, ratio
-real(r8), allocatable :: obs_temp(:), obs_qc(:), forward_temp(:)
+real(r8), allocatable :: obs_temp(:), forward_temp(:)
 real(r8)              :: obs_prior_mean, obs_prior_var, obs_val, obs_err_var
 real(r8)              :: rvalue(1)
-logical               :: do_outlier
+logical               :: do_outlier, good_forward_op
 
 
 ! Assume that mean and spread have been computed if needed???
@@ -1395,7 +1395,7 @@
 
 ! Compute mean and spread
 call compute_copy_mean_var(obs_ens_handle, &
-      1, ens_size, OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START)
+      1, ens_size, OBS_MEAN_START, OBS_VAR_START)
 
 
 ! Give the observation code a chance to alter the actual observation
@@ -1404,15 +1404,21 @@
 if (observations_updateable) then
   call update_observations_radar(obs_ens_handle, ens_size, seq, keys, prior_post, &
     obs_val_index, OBS_KEY_COPY, ens_mean_index, ens_spread_index, num_obs_in_set, &
-    OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
+    OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
     OBS_ERR_VAR_COPY, DART_qc_index, PRIOR_DIAG)
 endif
 
 
 ! At this point can compute outlier test and consolidate forward operator qc
 do j = 1, obs_ens_handle%my_num_vars
+   good_forward_op = .false.
+
+   ! find the min and max istatus values across all ensemble members.  these are
+   ! either set by dart code, or returned by the model-specific model_interpolate() 
+   ! routine, or by forward operator code in obs_def_xxx_mod files.
    forward_max = nint(maxval(forward_op_ens_handle%copies(1:ens_size, j)))
    forward_min = nint(minval(forward_op_ens_handle%copies(1:ens_size, j)))
+
    ! Now do a case statement to figure out what the qc result should be
    ! For prior, have to test for a bunch of stuff
    ! FIXME: note that this case statement doesn't cover every possibility;
@@ -1428,8 +1434,10 @@
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 4
       else if(forward_min == -1) then          ! Observation to be evaluated only
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 1
+         good_forward_op = .true.
       else if(forward_min == 0) then           ! All clear, assimilate this ob
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 0
+         good_forward_op = .true.
       ! FIXME: proposed enhancement - catchall for cases that we have not caught
       !else   ! 'should not happen'
       !   obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 9  ! inconsistent istatus codes
@@ -1443,8 +1451,8 @@
       ! only if it is still successful (assim or eval, 0 or 1), then check
       ! for failing outlier test.
       if(do_outlier .and. (obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) < 2)) then
-         obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START, j)
-         obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START, j)
+         obs_prior_mean = obs_ens_handle%copies(OBS_MEAN_START, j)
+         obs_prior_var = obs_ens_handle%copies(OBS_VAR_START, j)
          obs_val = obs_ens_handle%copies(OBS_VAL_COPY, j)
          obs_err_var = obs_ens_handle%copies(OBS_ERR_VAR_COPY, j)
          error = obs_prior_mean - obs_val
@@ -1474,7 +1482,26 @@
             obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 3
          endif
       endif
+      ! and for consistency, go through all the same tests as the prior, but
+      ! only use the results to set the good/bad forward op flag.
+      if ((forward_min == -99) .or. &       ! Failed prior qc in get_obs_ens
+          (forward_min == -2)  .or. &       ! Observation not used via namelist
+          (forward_max > 0)) then           ! At least one forward operator failed
+            continue; 
+      else if(forward_min == -1) then       ! Observation to be evaluated only
+         good_forward_op = .true.
+      else if(forward_min == 0) then        ! All clear, assimilate this ob
+         good_forward_op = .true.
+      endif
    endif
+
+   ! for either prior or posterior, if the forward operator failed,
+   ! reset the mean/var to missing_r8, regardless of the DART QC status
+   if (.not. good_forward_op) then
+      obs_ens_handle%copies(OBS_MEAN_START, j) = missing_r8
+      obs_ens_handle%copies(OBS_VAR_START,  j) = missing_r8
+   endif
+
 enddo
 
 ! PAR: NEED TO BE ABLE TO OUTPUT DETAILS OF FAILED FORWARD OBSEVATION OPERATORS
@@ -1485,60 +1512,32 @@
 call all_copies_to_all_vars(obs_ens_handle)
 
 ! allocate temp space for sending data
-allocate(obs_temp(num_obs_in_set), obs_qc(num_obs_in_set))
+allocate(obs_temp(num_obs_in_set))
 
-! set non-computed means & sd to missing values in the output diag file.   
-! for priors, the lowest good dart qc is 3, for posteriors its 1.
-! in either case, we want to see values that were only rejected because
-! of outlier threshold, so allow dart qc values 7 (and up?)
-if (prior_post == PRIOR_DIAG) then
-   low_qc_limit = 3
-else
-   low_qc_limit = 1
-endif
-high_qc_limit = 7
-
-! Update the qc global value
-! First get this copy on pe 0 and keep it around
-call get_copy(0, obs_ens_handle, OBS_GLOBAL_QC_COPY, obs_qc)
-! Only pe 0 gets to write the observations for this time
-if(my_task_id() == 0) then
-   ! Loop through the observations for this time
-   do j = 1, obs_ens_handle%num_vars
-      rvalue(1) = obs_qc(j)
-      call replace_qc(seq, keys(j), rvalue, DART_qc_index)
-   end do
-endif
-
 ! Update the ensemble mean
 ! Get this copy to process 0
-call get_copy(0, obs_ens_handle, OBS_PRIOR_MEAN_START, obs_temp)
+call get_copy(0, obs_ens_handle, OBS_MEAN_START, obs_temp)
 ! Only pe 0 gets to write the sequence
 if(my_task_id() == 0) then
      ! Loop through the observations for this time
      do j = 1, obs_ens_handle%num_vars
-      ! update the mean in each obs if it was assimilated ok.
-      if ((obs_qc(j) <= low_qc_limit) .or. (obs_qc(j) >= high_qc_limit)) then
-         rvalue(1) = obs_temp(j)
-      else
-         rvalue(1) = missing_r8    ! if qc bad, set to missing
-      endif
+      rvalue(1) = obs_temp(j)
       call replace_obs_values(seq, keys(j), rvalue, ens_mean_index)
      end do
   endif
 
 ! Update the ensemble spread
 ! Get this copy to process 0
-call get_copy(0, obs_ens_handle, OBS_PRIOR_VAR_START, obs_temp)
+call get_copy(0, obs_ens_handle, OBS_VAR_START, obs_temp)
 ! Only pe 0 gets to write the sequence
 if(my_task_id() == 0) then
    ! Loop through the observations for this time
    do j = 1, obs_ens_handle%num_vars
       ! update the spread in each obs
-      if ((obs_qc(j) <= low_qc_limit) .or. (obs_qc(j) >= high_qc_limit)) then
+      if (obs_temp(j) /= missing_r8) then
          rvalue(1) = sqrt(obs_temp(j))
       else
-         rvalue(1) = missing_r8    ! if qc bad, set to missing
+         rvalue(1) = obs_temp(j)
       endif
       call replace_obs_values(seq, keys(j), rvalue, ens_spread_index)
    end do
@@ -1563,8 +1562,19 @@
    endif
 end do
 
+! Update the qc global value
+call get_copy(0, obs_ens_handle, OBS_GLOBAL_QC_COPY, obs_temp)
+! Only pe 0 gets to write the observations for this time
+if(my_task_id() == 0) then
+   ! Loop through the observations for this time
+   do j = 1, obs_ens_handle%num_vars
+      rvalue(1) = obs_temp(j)
+      call replace_qc(seq, keys(j), rvalue, DART_qc_index)
+   end do
+endif
+
 ! clean up.
-deallocate(obs_temp, obs_qc)
+deallocate(obs_temp)
 
 end subroutine obs_space_diagnostics
 
@@ -1764,7 +1774,7 @@
 
 subroutine update_observations_radar(obs_ens_handle, ens_size, seq, keys, prior_post, &
       obs_val_index, OBS_KEY_COPY, ens_mean_index, ens_spread_index, num_obs_in_set, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
+      OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
       OBS_ERR_VAR_COPY, DART_qc_index, PRIOR_DIAG)
 
 
@@ -1786,7 +1796,7 @@
 integer,                 intent(in)    :: keys(num_obs_in_set), prior_post
 integer,                 intent(in)    :: ens_mean_index, ens_spread_index
 type(obs_sequence_type), intent(inout) :: seq
-integer,                 intent(in)    :: OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START
+integer,                 intent(in)    :: OBS_MEAN_START, OBS_VAR_START
 integer,                 intent(in)    :: OBS_GLOBAL_QC_COPY, OBS_VAL_COPY
 integer,                 intent(in)    :: OBS_ERR_VAR_COPY, DART_qc_index, PRIOR_DIAG
 
@@ -1829,7 +1839,7 @@
    call get_obs_def(observation, obs_def)
    obs_kind_ind = get_obs_kind(obs_def)
    if (obs_kind_ind == DOPPLER_RADIAL_VELOCITY) then
-      obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START, j)
+      obs_prior_mean = obs_ens_handle%copies(OBS_MEAN_START, j)
       obs_val = obs_ens_handle%copies(OBS_VAL_COPY, j)
       velkey = get_obs_def_key(obs_def)
       call get_obs_def_radial_vel(velkey, radarloc, beamdir, velnyquist)

Modified: DART/trunk/filter/filter.f90
===================================================================
--- DART/trunk/filter/filter.f90	2011-06-07 22:54:21 UTC (rev 4958)
+++ DART/trunk/filter/filter.f90	2011-06-08 14:55:02 UTC (rev 4959)
@@ -177,8 +177,8 @@
 integer                 :: ENS_MEAN_COPY, ENS_SD_COPY, PRIOR_INF_COPY, PRIOR_INF_SD_COPY
 integer                 :: POST_INF_COPY, POST_INF_SD_COPY
 integer                 :: OBS_VAL_COPY, OBS_ERR_VAR_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY
-integer                 :: OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END
-integer                 :: OBS_PRIOR_VAR_START, OBS_PRIOR_VAR_END, TOTAL_OBS_COPIES
+integer                 :: OBS_MEAN_START, OBS_MEAN_END
+integer                 :: OBS_VAR_START, OBS_VAR_END, TOTAL_OBS_COPIES
 integer                 :: input_qc_index, DART_qc_index
 integer                 :: mean_owner, mean_owners_index
 
@@ -241,10 +241,10 @@
 OBS_VAL_COPY         = ens_size + 2
 OBS_KEY_COPY         = ens_size + 3
 OBS_GLOBAL_QC_COPY   = ens_size + 4
-OBS_PRIOR_MEAN_START = ens_size + 5
-OBS_PRIOR_MEAN_END   = OBS_PRIOR_MEAN_START + num_groups - 1
-OBS_PRIOR_VAR_START  = OBS_PRIOR_MEAN_START + num_groups
-OBS_PRIOR_VAR_END    = OBS_PRIOR_VAR_START + num_groups - 1
+OBS_MEAN_START       = ens_size + 5
+OBS_MEAN_END         = OBS_MEAN_START + num_groups - 1
+OBS_VAR_START        = OBS_MEAN_START + num_groups
+OBS_VAR_END          = OBS_VAR_START + num_groups - 1
 
 ! Can't output more ensemble members than exist
 if(num_output_state_members > ens_size) num_output_state_members = ens_size
@@ -542,7 +542,7 @@
    call obs_space_diagnostics(obs_ens_handle, forward_op_ens_handle, ens_size, &
       seq, keys, PRIOR_DIAG, num_output_obs_members, in_obs_copy+1, &
       prior_obs_mean_index, prior_obs_spread_index, num_obs_in_set, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, &
+      OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
    call trace_message('After  observation space diagnostics')
   
@@ -559,8 +559,8 @@
       ens_size, num_groups, obs_val_index, prior_inflate, &
       ENS_MEAN_COPY, ENS_SD_COPY, &
       PRIOR_INF_COPY, PRIOR_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-      OBS_PRIOR_VAR_END, inflate_only = .false.)
+      OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+      OBS_VAR_END, inflate_only = .false.)
 
    call timestamp_message('After  observation assimilation')
    call     trace_message('After  observation assimilation')
@@ -577,8 +577,8 @@
       call smoother_assim(obs_ens_handle, seq, keys, ens_size, num_groups, &
          obs_val_index, ENS_MEAN_COPY, ENS_SD_COPY, &
          PRIOR_INF_COPY, PRIOR_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-         OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-         OBS_PRIOR_VAR_END)
+         OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+         OBS_VAR_END)
       call timestamp_message('After  smoother assimilation')
       call     trace_message('After  smoother assimilation')
    endif
@@ -650,7 +650,7 @@
    call obs_space_diagnostics(obs_ens_handle, forward_op_ens_handle, ens_size, &
       seq, keys, POSTERIOR_DIAG, num_output_obs_members, in_obs_copy+2, &
       posterior_obs_mean_index, posterior_obs_spread_index, num_obs_in_set, &
-      OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, &
+      OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
    
    call trace_message('After  posterior obs space diagnostics')
@@ -686,8 +686,8 @@
          call filter_assim(ens_handle, obs_ens_handle, seq, keys, ens_size, num_groups, &
             obs_val_index, post_inflate, ENS_MEAN_COPY, ENS_SD_COPY, &
             POST_INF_COPY, POST_INF_SD_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
-            OBS_PRIOR_MEAN_START, OBS_PRIOR_MEAN_END, OBS_PRIOR_VAR_START, &
-            OBS_PRIOR_VAR_END, inflate_only = .true.)
+            OBS_MEAN_START, OBS_MEAN_END, OBS_VAR_START, &
+            OBS_VAR_END, inflate_only = .true.)
 
          call all_copies_to_all_vars(ens_handle)
 
@@ -1308,7 +1308,7 @@
 subroutine obs_space_diagnostics(obs_ens_handle, forward_op_ens_handle, ens_size, &
    seq, keys, prior_post, num_output_members, members_index, &
    ens_mean_index, ens_spread_index, num_obs_in_set, &
-   OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
+   OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
    OBS_ERR_VAR_COPY, DART_qc_index)
 
 ! Do prior observation space diagnostics on the set of obs corresponding to keys
@@ -1320,17 +1320,17 @@
 integer,                 intent(in)    :: num_output_members, members_index
 integer,                 intent(in)    :: ens_mean_index, ens_spread_index
 type(obs_sequence_type), intent(inout) :: seq
-integer,                 intent(in)    :: OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START
+integer,                 intent(in)    :: OBS_MEAN_START, OBS_VAR_START
 integer,                 intent(in)    :: OBS_GLOBAL_QC_COPY, OBS_VAL_COPY
 integer,                 intent(in)    :: OBS_ERR_VAR_COPY, DART_qc_index
 
 integer               :: j, k, ens_offset, forward_min, forward_max
 integer               :: forward_unit, ivalue, low_qc_limit, high_qc_limit
 real(r8)              :: error, diff_sd, ratio
-real(r8), allocatable :: obs_temp(:), obs_qc(:), forward_temp(:)
+real(r8), allocatable :: obs_temp(:), forward_temp(:)
 real(r8)              :: obs_prior_mean, obs_prior_var, obs_val, obs_err_var
 real(r8)              :: rvalue(1)
-logical               :: do_outlier
+logical               :: do_outlier, good_forward_op
 
 
 ! Assume that mean and spread have been computed if needed???
@@ -1386,12 +1386,18 @@
 
 ! Compute mean and spread
 call compute_copy_mean_var(obs_ens_handle, &
-      1, ens_size, OBS_PRIOR_MEAN_START, OBS_PRIOR_VAR_START)
+      1, ens_size, OBS_MEAN_START, OBS_VAR_START)
 
 ! At this point can compute outlier test and consolidate forward operator qc
 do j = 1, obs_ens_handle%my_num_vars
+   good_forward_op = .false.
+
+   ! find the min and max istatus values across all ensemble members.  these are
+   ! either set by dart code, or returned by the model-specific model_interpolate() 
+   ! routine, or by forward operator code in obs_def_xxx_mod files.
    forward_max = nint(maxval(forward_op_ens_handle%copies(1:ens_size, j)))
    forward_min = nint(minval(forward_op_ens_handle%copies(1:ens_size, j)))
+
    ! Now do a case statement to figure out what the qc result should be
    ! For prior, have to test for a bunch of stuff
    ! FIXME: note that this case statement doesn't cover every possibility;
@@ -1407,8 +1413,10 @@
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 4
       else if(forward_min == -1) then          ! Observation to be evaluated only
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 1
+         good_forward_op = .true.
       else if(forward_min == 0) then           ! All clear, assimilate this ob
          obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 0
+         good_forward_op = .true.
       ! FIXME: proposed enhancement - catchall for cases that we have not caught
       !else   ! 'should not happen'
       !   obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 9  ! inconsistent istatus codes
@@ -1422,8 +1430,8 @@
       ! only if it is still successful (assim or eval, 0 or 1), then check
       ! for failing outlier test.
       if(do_outlier .and. (obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) < 2)) then
-         obs_prior_mean = obs_ens_handle%copies(OBS_PRIOR_MEAN_START, j)
-         obs_prior_var = obs_ens_handle%copies(OBS_PRIOR_VAR_START, j)
+         obs_prior_mean = obs_ens_handle%copies(OBS_MEAN_START, j)
+         obs_prior_var = obs_ens_handle%copies(OBS_VAR_START, j)
          obs_val = obs_ens_handle%copies(OBS_VAL_COPY, j)
          obs_err_var = obs_ens_handle%copies(OBS_ERR_VAR_COPY, j)
          error = obs_prior_mean - obs_val
@@ -1453,7 +1461,26 @@
             obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 3
          endif
       endif
+      ! and for consistency, go through all the same tests as the prior, but
+      ! only use the results to set the good/bad forward op flag.
+      if ((forward_min == -99) .or. &       ! Failed prior qc in get_obs_ens
+          (forward_min == -2)  .or. &       ! Observation not used via namelist
+          (forward_max > 0)) then           ! At least one forward operator failed
+            continue; 
+      else if(forward_min == -1) then       ! Observation to be evaluated only
+         good_forward_op = .true.
+      else if(forward_min == 0) then        ! All clear, assimilate this ob
+         good_forward_op = .true.
+      endif
    endif
+
+   ! for either prior or posterior, if the forward operator failed,
+   ! reset the mean/var to missing_r8, regardless of the DART QC status
+   if (.not. good_forward_op) then
+      obs_ens_handle%copies(OBS_MEAN_START, j) = missing_r8
+      obs_ens_handle%copies(OBS_VAR_START,  j) = missing_r8
+   endif
+
 enddo
 
 ! PAR: NEED TO BE ABLE TO OUTPUT DETAILS OF FAILED FORWARD OBSEVATION OPERATORS
@@ -1464,60 +1491,32 @@
 call all_copies_to_all_vars(obs_ens_handle)
 
 ! allocate temp space for sending data
-allocate(obs_temp(num_obs_in_set), obs_qc(num_obs_in_set))
+allocate(obs_temp(num_obs_in_set))
 
-! set non-computed means & sd to missing values in the output diag file.   
-! for priors, the lowest good dart qc is 3, for posteriors its 1.
-! in either case, we want to see values that were only rejected because
-! of outlier threshold, so allow dart qc values 7 (and up?)
-if (prior_post == PRIOR_DIAG) then
-   low_qc_limit = 3
-else
-   low_qc_limit = 1
-endif
-high_qc_limit = 7
-
-! Update the qc global value
-! First get this copy on pe 0 and keep it around
-call get_copy(0, obs_ens_handle, OBS_GLOBAL_QC_COPY, obs_qc)
-! Only pe 0 gets to write the observations for this time
-if(my_task_id() == 0) then
-   ! Loop through the observations for this time
-   do j = 1, obs_ens_handle%num_vars
-      rvalue(1) = obs_qc(j)
-      call replace_qc(seq, keys(j), rvalue, DART_qc_index)
-   end do
-endif
-
 ! Update the ensemble mean
 ! Get this copy to process 0
-call get_copy(0, obs_ens_handle, OBS_PRIOR_MEAN_START, obs_temp)
+call get_copy(0, obs_ens_handle, OBS_MEAN_START, obs_temp)
 ! Only pe 0 gets to write the sequence
 if(my_task_id() == 0) then
      ! Loop through the observations for this time
      do j = 1, obs_ens_handle%num_vars
-      ! update the mean in each obs if it was assimilated ok.
-      if ((obs_qc(j) <= low_qc_limit) .or. (obs_qc(j) >= high_qc_limit)) then
-         rvalue(1) = obs_temp(j)
-      else
-         rvalue(1) = missing_r8    ! if qc bad, set to missing
-      endif
+      rvalue(1) = obs_temp(j)
       call replace_obs_values(seq, keys(j), rvalue, ens_mean_index)
      end do
   endif
 
 ! Update the ensemble spread
 ! Get this copy to process 0
-call get_copy(0, obs_ens_handle, OBS_PRIOR_VAR_START, obs_temp)
+call get_copy(0, obs_ens_handle, OBS_VAR_START, obs_temp)
 ! Only pe 0 gets to write the sequence
 if(my_task_id() == 0) then
    ! Loop through the observations for this time
    do j = 1, obs_ens_handle%num_vars
       ! update the spread in each obs
-      if ((obs_qc(j) <= low_qc_limit) .or. (obs_qc(j) >= high_qc_limit)) then
+      if (obs_temp(j) /= missing_r8) then
          rvalue(1) = sqrt(obs_temp(j))
       else
-         rvalue(1) = missing_r8    ! if qc bad, set to missing
+         rvalue(1) = obs_temp(j)
       endif
       call replace_obs_values(seq, keys(j), rvalue, ens_spread_index)
    end do
@@ -1542,8 +1541,19 @@
    endif
 end do
 
+! Update the qc global value
+call get_copy(0, obs_ens_handle, OBS_GLOBAL_QC_COPY, obs_temp)
+! Only pe 0 gets to write the observations for this time
+if(my_task_id() == 0) then
+   ! Loop through the observations for this time
+   do j = 1, obs_ens_handle%num_vars
+      rvalue(1) = obs_temp(j)
+      call replace_qc(seq, keys(j), rvalue, DART_qc_index)
+   end do
+endif
+
 ! clean up.
-deallocate(obs_temp, obs_qc)
+deallocate(obs_temp)
 
 end subroutine obs_space_diagnostics
 


More information about the Dart-dev mailing list