[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