[Dart-dev] [6011] DART/branches/development/filter: Moved the test for the outlier_threshold into a separate

nancy at ucar.edu nancy at ucar.edu
Thu Mar 21 11:49:47 MDT 2013


Revision: 6011
Author:   nancy
Date:     2013-03-21 11:49:46 -0600 (Thu, 21 Mar 2013)
Log Message:
-----------
Moved the test for the outlier_threshold into a separate
subroutine at the end of the file.  Added comments about
how users can customize this subroutine if they want a
different threshold based on location or observation type.
There is a new namelist item, enable_special_outlier_code,
which can be turned on and off at run time to make it easy
to run comparisons between the original simple threshold
and any more complex handling added by the user.

Also added a few new lines to the log file - if inflation
damping is on, echo that value.  Echo when calling the
initial and end routines in the model_mod, to help isolate
problems if the code fails.  Also call end_model() for
the first time. It should have been called all along
and never was.  our bad.  Make the qc threshold code
a bit more modular by not referencing a global var.

Identical changes to filter.dopplerfold.f90 to keep
in sync with filter; changes to the default nml and
html pages to document the new namelist item.

All this being said, there should be no changes to
the filter results from this update; there will be
a few new lines in the log files but that's all.

Modified Paths:
--------------
    DART/branches/development/filter/filter.dopplerfold.f90
    DART/branches/development/filter/filter.f90
    DART/branches/development/filter/filter.html
    DART/branches/development/filter/filter.nml

-------------- next part --------------
Modified: DART/branches/development/filter/filter.dopplerfold.f90
===================================================================
--- DART/branches/development/filter/filter.dopplerfold.f90	2013-03-15 19:31:17 UTC (rev 6010)
+++ DART/branches/development/filter/filter.dopplerfold.f90	2013-03-21 17:49:46 UTC (rev 6011)
@@ -31,7 +31,7 @@
                                  open_file, close_file, do_nml_file, do_nml_term
 use assim_model_mod,      only : static_init_assim_model, get_model_size,                    &
                                  netcdf_file_type, init_diag_output, finalize_diag_output,   & 
-                                 aoutput_diagnostics, ens_mean_for_model
+                                 aoutput_diagnostics, ens_mean_for_model, end_assim_model
 use assim_tools_mod,      only : filter_assim, set_assim_tools_trace
 use obs_model_mod,        only : move_ahead, advance_state, set_obs_model_trace
 use ensemble_manager_mod, only : init_ensemble_manager, end_ensemble_manager,                &
@@ -101,6 +101,7 @@
 integer  :: output_interval     = 1
 integer  :: num_groups          = 1
 real(r8) :: outlier_threshold   = -1.0_r8
+logical  :: enable_special_outlier_code = .false.
 real(r8) :: input_qc_threshold  = 3.0_r8
 logical  :: output_forward_op_errors = .false.
 logical  :: output_timestamps        = .false.
@@ -137,7 +138,7 @@
    start_from_restart, output_restart, obs_sequence_in_name, obs_sequence_out_name, &
    restart_in_file_name, restart_out_file_name, init_time_days, init_time_seconds,  &
    first_obs_days, first_obs_seconds, last_obs_days, last_obs_seconds,              &
-   obs_window_days, obs_window_seconds,                                             &
+   obs_window_days, obs_window_seconds, enable_special_outlier_code,                &
    num_output_state_members, num_output_obs_members, output_restart_mean,           &
    output_interval, num_groups, outlier_threshold, trace_execution,                 &
    input_qc_threshold, output_forward_op_errors, output_timestamps,                 &
@@ -215,6 +216,10 @@
    call error_handler(E_ERR,'filter_main', msgstring, source, revision, revdate)
 endif
 
+! informational message to log
+write(msgstring, '(A,I5)') 'running with an ensemble size of ', ens_size
+call error_handler(E_MSG,'filter:', msgstring, source, revision, revdate)
+
 ! See if smoothing is turned on
 ds = do_smoothing()
 
@@ -230,6 +235,12 @@
    endif
 end do
 
+! if doing something special with outlier threshold, say so
+if (enable_special_outlier_code) then
+   call error_handler(E_MSG,'filter:', 'special outlier threshold handling enabled', &
+      source, revision, revdate)
+endif
+
 ! Observation space inflation for posterior not currently supported
 if(inf_flavor(2) == 1) call error_handler(E_ERR, 'filter_main', &
    'Posterior observation space inflation (type 1) not supported', source, revision, revdate)
@@ -308,6 +319,17 @@
    inf_sd_initial(2), inf_lower_bound(2), inf_upper_bound(2), inf_sd_lower_bound(2),  &
    ens_handle, POST_INF_COPY, POST_INF_SD_COPY, 'Posterior')
 
+if (do_output()) then
+   if (inf_flavor(1) > 0 .and. inf_damping(1) < 1.0_r8) then
+      write(msgstring, '(A,F12.6,A)') 'Prior inflation damping of ', inf_damping(1), ' will be used'
+      call error_handler(E_MSG,'filter:', msgstring)
+   endif
+   if (inf_flavor(2) > 0 .and. inf_damping(2) < 1.0_r8) then
+      write(msgstring, '(A,F12.6,A)') 'Posterior inflation damping of ', inf_damping(2), ' will be used'
+      call error_handler(E_MSG,'filter:', msgstring)
+   endif
+endif
+
 call trace_message('After  initializing inflation')
 
 call     trace_message('Before initializing output files')
@@ -421,6 +443,7 @@
       exit AdvanceTime
    endif
 
+
    ! if model state data not at required time, advance model
    if (curr_ens_time /= next_ens_time) then
       ! Advance the lagged distribution, if needed.
@@ -538,11 +561,13 @@
    if ((output_interval > 0) .and. &
        (time_step_number / output_interval * output_interval == time_step_number)) then
       call trace_message('Before prior state space diagnostics')
+      call timestamp_message('Before prior state space diagnostics')
       call filter_state_space_diagnostics(PriorStateUnit, ens_handle, &
          model_size, num_output_state_members, &
          output_state_mean_index, output_state_spread_index, &
          output_inflation, ens_mean, ENS_MEAN_COPY, ENS_SD_COPY, &
          prior_inflate, PRIOR_INF_COPY, PRIOR_INF_SD_COPY)
+      call timestamp_message('After  prior state space diagnostics')
       call trace_message('After  prior state space diagnostics')
    endif
   
@@ -641,6 +666,7 @@
    if ((output_interval > 0) .and. &
        (time_step_number / output_interval * output_interval == time_step_number)) then
       call trace_message('Before posterior state space diagnostics')
+      call timestamp_message('Before posterior state space diagnostics')
       call filter_state_space_diagnostics(PosteriorStateUnit, ens_handle, &
          model_size, num_output_state_members, output_state_mean_index, &
          output_state_spread_index, &
@@ -651,6 +677,7 @@
       call smoother_ss_diagnostics(model_size, num_output_state_members, &
          output_inflation, ens_mean, ENS_MEAN_COPY, ENS_SD_COPY, &
          POST_INF_COPY, POST_INF_SD_COPY)
+      call timestamp_message('After  posterior state space diagnostics')
       call trace_message('After  posterior state space diagnostics')
    endif
 
@@ -760,6 +787,11 @@
 if(ds) call smoother_write_restart(1, ens_size)
 call trace_message('After  writing state restart files if requested')
 
+! Give the model_mod code a chance to clean up. 
+call trace_message('Before end_model call')
+call end_assim_model()
+call trace_message('After  end_model call')
+
 call trace_message('Before ensemble and obs memory cleanup')
 call end_ensemble_manager(ens_handle)
 
@@ -919,7 +951,9 @@
 call static_init_obs_sequence()
 
 ! Initialize the model class data now that obs_sequence is all set up
+call trace_message('Before init_model call')
 call static_init_assim_model()
+call trace_message('After  init_model call')
 
 end subroutine filter_initialize_modules_used
 
@@ -1238,7 +1272,7 @@
    ! If it is bad, set forward operator status value to -99 and return missing_r8 for obs_value
 
    ! PAR THIS SUBROUTINE SHOULD EVENTUALLY GO IN THE QUALITY CONTROL MODULE
-   if(.not. input_qc_ok(input_qc(1))) then
+   if(.not. input_qc_ok(input_qc(1), input_qc_threshold)) then
       ! The forward operator value is set to -99 if prior qc was failed
       forward_op_ens_handle%vars(j, :) = -99
       do k=1, my_num_copies
@@ -1250,9 +1284,6 @@
             obs_ens_handle%vars(j, k) = missing_r8
         endif
       enddo
-      ! previous code was at various times one of these equivalent lines:
-      !obs_ens_handle%vars(j, 1:my_num_copies) = missing_r8
-      !obs_ens_handle%vars(j, :) = missing_r8 
       ! No need to do anything else for a failed observation
       cycle ALL_OBSERVATIONS
    endif
@@ -1266,14 +1297,7 @@
       global_ens_index = obs_ens_handle%my_copies(k)
       ! If I have a copy that is a standard ensemble member, compute expected value
       if(global_ens_index <= ens_size) then
-         ! Compute the expected observation value; direct access to storage
-         ! Debug: The PGI 6.1.3 compiler was making internal copies of the array sections
-         ! very slowly, causing the filter process to take > 12 hours (vs minutes).
-         ! Try making local copies/non-sections to convince the compiler to run fast.
-         ! The original code was:
-         !    call get_expected_obs(seq, keys(j:j), ens_handle%vars(:, k), &
-         !       obs_ens_handle%vars(j:j, k), istatus, assimilate_this_ob, evaluate_this_ob)
-         ! and keys is intent in only, vars intent out only.
+         ! temporaries to avoid passing array sections which was slow on PGI compiler
          thiskey(1) = keys(j)
          call get_expected_obs(seq, thiskey, &
             global_ens_index, ens_handle%vars(:, k), ens_handle%time(1), &
@@ -1344,7 +1368,7 @@
 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, good_forward_op
+logical               :: do_outlier, good_forward_op, failed
 
 
 ! Assume that mean and spread have been computed if needed???
@@ -1462,16 +1486,24 @@
          obs_err_var = obs_ens_handle%copies(OBS_ERR_VAR_COPY, j)
          error = obs_prior_mean - obs_val
          diff_sd = sqrt(obs_prior_var + obs_err_var)
-         ratio = abs(error / diff_sd)
-         if(ratio > outlier_threshold) then
-            ! FIXME: proposed enhancement to dart qc values - differentiate
-            ! between obs thrown out for being outlier - assim vs eval.
-            !if(obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) == 0) then 
-               obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 7    ! would have been assim
-            !else
-            !   obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 8    ! was eval only
-            !endif
+         if (diff_sd /= 0.0_r8) then
+            ratio = abs(error / diff_sd)
+         else
+            ratio = 0.0_r8
          endif
+
+         ! if special handling requested, pass in the outlier ratio for this obs,
+         ! the default outlier threshold value, and enough info to extract the specific 
+         ! obs type for this obs.
+         ! the function should return .true. if this is an outlier, .false. if it is ok.
+         if (enable_special_outlier_code) then
+            failed = failed_outlier(ratio, outlier_threshold, obs_ens_handle, &
+                                    OBS_KEY_COPY, j, seq)
+         else 
+            failed = (ratio > outlier_threshold)
+         endif
+
+         if (failed) obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 7 
       endif
 
    else
@@ -1585,10 +1617,10 @@
 
 !-------------------------------------------------------------------------
 
-function input_qc_ok(input_qc)
+function input_qc_ok(input_qc, threshold)
 
 logical              :: input_qc_ok
-real(r8), intent(in) :: input_qc
+real(r8), intent(in) :: input_qc, threshold
 
 ! Do checks on input_qc value with namelist control
 ! Should eventually go in qc module
@@ -1599,10 +1631,7 @@
 ! e.g. to keep only obs with a data qc of 0, set the 
 ! threshold to 0.
 
-! To exclude negative qc values, comment in the following line instead
-! of the existing code.
-! if((nint(input_qc) <= nint(input_qc_threshold)) .and. (nint(input_qc) >= 0)) then
-if(nint(input_qc) <= nint(input_qc_threshold)) then
+if(nint(input_qc) <= nint(threshold)) then
    input_qc_ok = .true.
 else
    input_qc_ok = .false.
@@ -1773,7 +1802,83 @@
 
 !-------------------------------------------------------------------------
 
+function failed_outlier(ratio, outlier_threshold, obs_ens_handle, OBS_KEY_COPY, j, seq)
 
+! return true if the observation value is too far away from the ensemble mean
+! and should be rejected and not assimilated.
+
+use obs_def_mod, only : get_obs_kind
+use obs_kind_mod         ! this allows you to use all the types available
+
+
+real(r8),                intent(in) :: ratio
+real(r8),                intent(in) :: outlier_threshold
+type(ensemble_type),     intent(in) :: obs_ens_handle
+integer,                 intent(in) :: OBS_KEY_COPY
+integer,                 intent(in) :: j
+type(obs_sequence_type), intent(in) :: seq
+logical                             :: failed_outlier
+
+! the default test is:  if (ratio > outlier_threshold) failed_outlier = .true.
+! but you can add code here to do different tests for different observation 
+! types.  this function is only called if this namelist item is set to true:
+!  enable_special_outlier_code = .true.
+! in the &filter_nml namelist.  it is intended to be customized by the user.
+
+integer :: this_obs_key, this_obs_type
+type(obs_def_type) :: obs_def
+type(obs_type) :: observation
+logical :: first_time = .true.
+
+! make sure there's space to hold the observation.   this is a memory
+! leak in that we never release this space, but we only allocate it
+! one time so it doesn't grow.  if you are going to access the values
+! of the observation or the qc values, you must change the 0, 0 below
+! to match what's in your obs_seq file.  the example code below uses
+! only things in the obs_def derived type and so doesn't need space
+! allocated for either copies or qcs.
+if (first_time) then
+   call init_obs(observation, 0, 0)
+   first_time = .false.
+endif
+
+! if you want to do something different based on the observation specific type:
+
+this_obs_key = obs_ens_handle%copies(OBS_KEY_COPY, j)
+call get_obs_from_key(seq, this_obs_key, observation)
+call get_obs_def(observation, obs_def)
+this_obs_type = get_obs_kind(obs_def)
+
+! note that this example uses the specific type (e.g. RADIOSONDE_TEMPERATURE)
+! to make decisions.  you have the observation so any other part (e.g. the
+! time, the value, the error) is available to you as well.
+
+select case (this_obs_type)
+
+! example of specifying a different threshold value for one obs type:
+!   case (RADIOSONDE_TEMPERATURE)
+!      if (ratio > some_other_value) then
+!         failed_outlier = .true.
+!      else
+!         failed_outlier = .false.
+!      endif
+
+! accept all values of this observation type no matter how far
+! from the ensemble mean:
+!   case (AIRCRAFT_U_WIND_COMPONENT, AIRCRAFT_V_WIND_COMPONENT)
+!      failed_outlier = .false.
+
+   case default
+      if (ratio > outlier_threshold) then
+         failed_outlier = .true.
+      else
+         failed_outlier = .false.
+      endif
+
+end select
+
+end function failed_outlier
+
 !-------------------------------------------------------------------------
 !-------------------------------------------------------------------------
 ! special for doppler radar unfolding

Modified: DART/branches/development/filter/filter.f90
===================================================================
--- DART/branches/development/filter/filter.f90	2013-03-15 19:31:17 UTC (rev 6010)
+++ DART/branches/development/filter/filter.f90	2013-03-21 17:49:46 UTC (rev 6011)
@@ -31,7 +31,7 @@
                                  open_file, close_file, do_nml_file, do_nml_term
 use assim_model_mod,      only : static_init_assim_model, get_model_size,                    &
                                  netcdf_file_type, init_diag_output, finalize_diag_output,   & 
-                                 aoutput_diagnostics, ens_mean_for_model
+                                 aoutput_diagnostics, ens_mean_for_model, end_assim_model
 use assim_tools_mod,      only : filter_assim, set_assim_tools_trace
 use obs_model_mod,        only : move_ahead, advance_state, set_obs_model_trace
 use ensemble_manager_mod, only : init_ensemble_manager, end_ensemble_manager,                &
@@ -101,6 +101,7 @@
 integer  :: output_interval     = 1
 integer  :: num_groups          = 1
 real(r8) :: outlier_threshold   = -1.0_r8
+logical  :: enable_special_outlier_code = .false.
 real(r8) :: input_qc_threshold  = 3.0_r8
 logical  :: output_forward_op_errors = .false.
 logical  :: output_timestamps        = .false.
@@ -137,7 +138,7 @@
    start_from_restart, output_restart, obs_sequence_in_name, obs_sequence_out_name, &
    restart_in_file_name, restart_out_file_name, init_time_days, init_time_seconds,  &
    first_obs_days, first_obs_seconds, last_obs_days, last_obs_seconds,              &
-   obs_window_days, obs_window_seconds,                                             &
+   obs_window_days, obs_window_seconds, enable_special_outlier_code,                &
    num_output_state_members, num_output_obs_members, output_restart_mean,           &
    output_interval, num_groups, outlier_threshold, trace_execution,                 &
    input_qc_threshold, output_forward_op_errors, output_timestamps,                 &
@@ -212,8 +213,8 @@
 endif
 
 ! informational message to log
-write(msgstring, *) 'running with an ensemble size of ', ens_size
-call error_handler(E_MSG,'filter:', msgstring)
+write(msgstring, '(A,I5)') 'running with an ensemble size of ', ens_size
+call error_handler(E_MSG,'filter:', msgstring, source, revision, revdate)
 
 ! See if smoothing is turned on
 ds = do_smoothing()
@@ -230,6 +231,12 @@
    endif
 end do
 
+! if doing something special with outlier threshold, say so
+if (enable_special_outlier_code) then
+   call error_handler(E_MSG,'filter:', 'special outlier threshold handling enabled', &
+      source, revision, revdate)
+endif
+
 ! Observation space inflation for posterior not currently supported
 if(inf_flavor(2) == 1) call error_handler(E_ERR, 'filter_main', &
    'Posterior observation space inflation (type 1) not supported', source, revision, revdate)
@@ -308,6 +315,17 @@
    inf_sd_initial(2), inf_lower_bound(2), inf_upper_bound(2), inf_sd_lower_bound(2),  &
    ens_handle, POST_INF_COPY, POST_INF_SD_COPY, 'Posterior')
 
+if (do_output()) then
+   if (inf_flavor(1) > 0 .and. inf_damping(1) < 1.0_r8) then
+      write(msgstring, '(A,F12.6,A)') 'Prior inflation damping of ', inf_damping(1), ' will be used'
+      call error_handler(E_MSG,'filter:', msgstring)
+   endif
+   if (inf_flavor(2) > 0 .and. inf_damping(2) < 1.0_r8) then
+      write(msgstring, '(A,F12.6,A)') 'Posterior inflation damping of ', inf_damping(2), ' will be used'
+      call error_handler(E_MSG,'filter:', msgstring)
+   endif
+endif
+
 call trace_message('After  initializing inflation')
 
 call     trace_message('Before initializing output files')
@@ -421,6 +439,7 @@
       exit AdvanceTime
    endif
 
+
    ! if model state data not at required time, advance model
    if (curr_ens_time /= next_ens_time) then
       ! Advance the lagged distribution, if needed.
@@ -538,11 +557,13 @@
    if ((output_interval > 0) .and. &
        (time_step_number / output_interval * output_interval == time_step_number)) then
       call trace_message('Before prior state space diagnostics')
+      call timestamp_message('Before prior state space diagnostics')
       call filter_state_space_diagnostics(PriorStateUnit, ens_handle, &
          model_size, num_output_state_members, &
          output_state_mean_index, output_state_spread_index, &
          output_inflation, ens_mean, ENS_MEAN_COPY, ENS_SD_COPY, &
          prior_inflate, PRIOR_INF_COPY, PRIOR_INF_SD_COPY)
+      call timestamp_message('After  prior state space diagnostics')
       call trace_message('After  prior state space diagnostics')
    endif
   
@@ -550,6 +571,7 @@
    ! Do prior observation space diagnostics and associated quality control
    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, &
+      OBS_KEY_COPY, &                                 ! new
       prior_obs_mean_index, prior_obs_spread_index, num_obs_in_set, &
       OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
@@ -640,6 +662,7 @@
    if ((output_interval > 0) .and. &
        (time_step_number / output_interval * output_interval == time_step_number)) then
       call trace_message('Before posterior state space diagnostics')
+      call timestamp_message('Before posterior state space diagnostics')
       call filter_state_space_diagnostics(PosteriorStateUnit, ens_handle, &
          model_size, num_output_state_members, output_state_mean_index, &
          output_state_spread_index, &
@@ -650,6 +673,7 @@
       call smoother_ss_diagnostics(model_size, num_output_state_members, &
          output_inflation, ens_mean, ENS_MEAN_COPY, ENS_SD_COPY, &
          POST_INF_COPY, POST_INF_SD_COPY)
+      call timestamp_message('After  posterior state space diagnostics')
       call trace_message('After  posterior state space diagnostics')
    endif
 
@@ -658,6 +682,7 @@
    ! Do posterior observation space diagnostics
    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, &
+      OBS_KEY_COPY, &                             ! new
       posterior_obs_mean_index, posterior_obs_spread_index, num_obs_in_set, &
       OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, &
       OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
@@ -758,6 +783,11 @@
 if(ds) call smoother_write_restart(1, ens_size)
 call trace_message('After  writing state restart files if requested')
 
+! Give the model_mod code a chance to clean up. 
+call trace_message('Before end_model call')
+call end_assim_model()
+call trace_message('After  end_model call')
+
 call trace_message('Before ensemble and obs memory cleanup')
 call end_ensemble_manager(ens_handle)
 
@@ -917,7 +947,9 @@
 call static_init_obs_sequence()
 
 ! Initialize the model class data now that obs_sequence is all set up
+call trace_message('Before init_model call')
 call static_init_assim_model()
+call trace_message('After  init_model call')
 
 end subroutine filter_initialize_modules_used
 
@@ -1236,7 +1268,7 @@
    ! If it is bad, set forward operator status value to -99 and return missing_r8 for obs_value
 
    ! PAR THIS SUBROUTINE SHOULD EVENTUALLY GO IN THE QUALITY CONTROL MODULE
-   if(.not. input_qc_ok(input_qc(1))) then
+   if(.not. input_qc_ok(input_qc(1), input_qc_threshold)) then
       ! The forward operator value is set to -99 if prior qc was failed
       forward_op_ens_handle%vars(j, :) = -99
       do k=1, my_num_copies
@@ -1248,9 +1280,6 @@
             obs_ens_handle%vars(j, k) = missing_r8
         endif
       enddo
-      ! previous code was at various times one of these equivalent lines:
-      !obs_ens_handle%vars(j, 1:my_num_copies) = missing_r8
-      !obs_ens_handle%vars(j, :) = missing_r8 
       ! No need to do anything else for a failed observation
       cycle ALL_OBSERVATIONS
    endif
@@ -1264,14 +1293,7 @@
       global_ens_index = obs_ens_handle%my_copies(k)
       ! If I have a copy that is a standard ensemble member, compute expected value
       if(global_ens_index <= ens_size) then
-         ! Compute the expected observation value; direct access to storage
-         ! Debug: The PGI 6.1.3 compiler was making internal copies of the array sections
-         ! very slowly, causing the filter process to take > 12 hours (vs minutes).
-         ! Try making local copies/non-sections to convince the compiler to run fast.
-         ! The original code was:
-         !    call get_expected_obs(seq, keys(j:j), ens_handle%vars(:, k), &
-         !       obs_ens_handle%vars(j:j, k), istatus, assimilate_this_ob, evaluate_this_ob)
-         ! and keys is intent in only, vars intent out only.
+         ! temporaries to avoid passing array sections which was slow on PGI compiler
          thiskey(1) = keys(j)
          call get_expected_obs(seq, thiskey, &
             global_ens_index, ens_handle%vars(:, k), ens_handle%time(1), &
@@ -1316,6 +1338,7 @@
 
 subroutine obs_space_diagnostics(obs_ens_handle, forward_op_ens_handle, ens_size, &
    seq, keys, prior_post, num_output_members, members_index, &
+   OBS_KEY_COPY, &
    ens_mean_index, ens_spread_index, num_obs_in_set, &
    OBS_MEAN_START, OBS_VAR_START, OBS_GLOBAL_QC_COPY, OBS_VAL_COPY, &
    OBS_ERR_VAR_COPY, DART_qc_index)
@@ -1327,6 +1350,7 @@
 integer,                 intent(in)    :: num_obs_in_set
 integer,                 intent(in)    :: keys(num_obs_in_set), prior_post
 integer,                 intent(in)    :: num_output_members, members_index
+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_MEAN_START, OBS_VAR_START
@@ -1339,7 +1363,7 @@
 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, good_forward_op
+logical               :: do_outlier, good_forward_op, failed
 
 
 ! Assume that mean and spread have been computed if needed???
@@ -1445,16 +1469,24 @@
          obs_err_var = obs_ens_handle%copies(OBS_ERR_VAR_COPY, j)
          error = obs_prior_mean - obs_val
          diff_sd = sqrt(obs_prior_var + obs_err_var)
-         ratio = abs(error / diff_sd)
-         if(ratio > outlier_threshold) then
-            ! FIXME: proposed enhancement to dart qc values - differentiate
-            ! between obs thrown out for being outlier - assim vs eval.
-            !if(obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) == 0) then 
-               obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 7    ! would have been assim
-            !else
-            !   obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 8    ! was eval only
-            !endif
+         if (diff_sd /= 0.0_r8) then
+            ratio = abs(error / diff_sd)
+         else
+            ratio = 0.0_r8
          endif
+
+         ! if special handling requested, pass in the outlier ratio for this obs,
+         ! the default outlier threshold value, and enough info to extract the specific 
+         ! obs type for this obs.
+         ! the function should return .true. if this is an outlier, .false. if it is ok.
+         if (enable_special_outlier_code) then
+            failed = failed_outlier(ratio, outlier_threshold, obs_ens_handle, &
+                                    OBS_KEY_COPY, j, seq)
+         else 
+            failed = (ratio > outlier_threshold)
+         endif
+
+         if (failed) obs_ens_handle%copies(OBS_GLOBAL_QC_COPY, j) = 7 
       endif
 
    else
@@ -1568,10 +1600,10 @@
 
 !-------------------------------------------------------------------------
 
-function input_qc_ok(input_qc)
+function input_qc_ok(input_qc, threshold)
 
 logical              :: input_qc_ok
-real(r8), intent(in) :: input_qc
+real(r8), intent(in) :: input_qc, threshold
 
 ! Do checks on input_qc value with namelist control
 ! Should eventually go in qc module
@@ -1582,10 +1614,7 @@
 ! e.g. to keep only obs with a data qc of 0, set the 
 ! threshold to 0.
 
-! To exclude negative qc values, comment in the following line instead
-! of the existing code.
-! if((nint(input_qc) <= nint(input_qc_threshold)) .and. (nint(input_qc) >= 0)) then
-if(nint(input_qc) <= nint(input_qc_threshold)) then
+if(nint(input_qc) <= nint(threshold)) then
    input_qc_ok = .true.
 else
    input_qc_ok = .false.
@@ -1756,4 +1785,83 @@
 
 !-------------------------------------------------------------------------
 
+function failed_outlier(ratio, outlier_threshold, obs_ens_handle, OBS_KEY_COPY, j, seq)
+
+! return true if the observation value is too far away from the ensemble mean
+! and should be rejected and not assimilated.
+
+use obs_def_mod, only : get_obs_kind
+use obs_kind_mod         ! this allows you to use all the types available
+
+
+real(r8),                intent(in) :: ratio
+real(r8),                intent(in) :: outlier_threshold
+type(ensemble_type),     intent(in) :: obs_ens_handle
+integer,                 intent(in) :: OBS_KEY_COPY
+integer,                 intent(in) :: j
+type(obs_sequence_type), intent(in) :: seq
+logical                             :: failed_outlier
+
+! the default test is:  if (ratio > outlier_threshold) failed_outlier = .true.
+! but you can add code here to do different tests for different observation 
+! types.  this function is only called if this namelist item is set to true:
+!  enable_special_outlier_code = .true.
+! in the &filter_nml namelist.  it is intended to be customized by the user.
+
+integer :: this_obs_key, this_obs_type
+type(obs_def_type) :: obs_def
+type(obs_type) :: observation
+logical :: first_time = .true.
+
+! make sure there's space to hold the observation.   this is a memory
+! leak in that we never release this space, but we only allocate it
+! one time so it doesn't grow.  if you are going to access the values
+! of the observation or the qc values, you must change the 0, 0 below
+! to match what's in your obs_seq file.  the example code below uses
+! only things in the obs_def derived type and so doesn't need space
+! allocated for either copies or qcs.
+if (first_time) then
+   call init_obs(observation, 0, 0)
+   first_time = .false.
+endif
+
+! if you want to do something different based on the observation specific type:
+
+this_obs_key = obs_ens_handle%copies(OBS_KEY_COPY, j)
+call get_obs_from_key(seq, this_obs_key, observation)
+call get_obs_def(observation, obs_def)
+this_obs_type = get_obs_kind(obs_def)
+
+! note that this example uses the specific type (e.g. RADIOSONDE_TEMPERATURE)
+! to make decisions.  you have the observation so any other part (e.g. the
+! time, the value, the error) is available to you as well.
+
+select case (this_obs_type)
+
+! example of specifying a different threshold value for one obs type:
+!   case (RADIOSONDE_TEMPERATURE)
+!      if (ratio > some_other_value) then
+!         failed_outlier = .true.
+!      else
+!         failed_outlier = .false.
+!      endif
+
+! accept all values of this observation type no matter how far
+! from the ensemble mean:
+!   case (AIRCRAFT_U_WIND_COMPONENT, AIRCRAFT_V_WIND_COMPONENT)
+!      failed_outlier = .false.
+
+   case default
+      if (ratio > outlier_threshold) then
+         failed_outlier = .true.
+      else
+         failed_outlier = .false.
+      endif
+
+end select
+
+end function failed_outlier
+
+!-------------------------------------------------------------------------
+
 end program filter

Modified: DART/branches/development/filter/filter.html
===================================================================
--- DART/branches/development/filter/filter.html	2013-03-15 19:31:17 UTC (rev 6010)
+++ DART/branches/development/filter/filter.html	2013-03-21 17:49:46 UTC (rev 6011)
@@ -642,6 +642,7 @@
   num_groups,
   input_qc_threshold,
   outlier_threshold,
+  enable_special_outlier_code,
   output_forward_op_errors,
   output_restart_mean,
   output_timestamps,
@@ -817,6 +818,21 @@
                        sd's from observation. Negative means no check. 
                        Default: -1.0_r8</TD></TR>
 
+<TR><!--contents--><TD valign=top>enable_special_outlier_code</TD>
+    <!--  type  --><TD valign=top>logical</TD>
+    <!--descript--><TD>If true call a subroutine which can be customized by the
+                       user to do more elaborate outlier thresholding tests.
+                       See <em class=code>failed_outlier()</em> near the bottom
+                       of <em class=file>filter.f90</em> for where to add the
+                       custom code, and for commented out examples of possible tests.
+                       Turning this flag on and off allows comparisons to be 
+                       made between the default outlier threshold code and any
+                       custom settings without having to recompile the code. 
+                       To change the outlier behavior code has to be changed
+                       in filter.f90.  As distributed, turning this flag on
+                       and off will make no difference in the results.
+                       Default: .false.</TD></TR>
+
 <TR><!--contents--><TD valign=top>input_qc_threshold</TD>
     <!--  type  --><TD valign=top>real(r8)</TD>
     <!--descript--><TD>Reject observation if incoming QC value exceeds

Modified: DART/branches/development/filter/filter.nml
===================================================================
--- DART/branches/development/filter/filter.nml	2013-03-15 19:31:17 UTC (rev 6010)
+++ DART/branches/development/filter/filter.nml	2013-03-21 17:49:46 UTC (rev 6011)
@@ -20,6 +20,7 @@
    num_groups               = 1,
    input_qc_threshold       =  3.0,
    outlier_threshold        = -1.0,
+   enable_special_outlier_code = .false.,
    output_forward_op_errors = .false.,
    output_restart_mean      = .false.,
    output_timestamps        = .false.,


More information about the Dart-dev mailing list