[Dart-dev] [6200] DART/branches/development:
nancy at ucar.edu
nancy at ucar.edu
Thu May 30 12:47:40 MDT 2013
Revision: 6200
Author: nancy
Date: 2013-05-30 12:47:39 -0600 (Thu, 30 May 2013)
Log Message:
-----------
fixes for the serious specific types/generic kinds mixup.
now in all get_close_obs() calls the base is always passed
as a specific type, and the target list is always a list
of generic kinds. where possible i changed the names of
things to be accurate.
in assim_tools there is new code to check for either
state vector items with 0 spread, or observations with
0 error variance. these are handled before the filter_kind
routines are called, so 1) they are all handled in a consistent
way and 2) the specific filter_kind routines don't have to do
error checking (some used to, some didn't - one user reported
a NaN error in filter kind 3 when he had 0 spread in the
ensemble members). these cases will now all be caught and
handled consistently.
added a new namelist item to allow disabling of the buffering
for searches when observations in the same location. for
things like U, V winds, or T, U, V readings from radiosondes
at mandatory levels where several successive observations are
at identical locations, the buffering can be a large time savings.
but if the model_mod has get_close_obs() code which depends
on looking at the specific type of the current obs, you can't
buffer previous searches at the same location. the buffering
is on by default, but the user can disable it now.
extensive internal changes in the DEFAULT_obs_kind_mod.F90
to try to separate out the specific types and generic kinds
inside the routines. next release we may be forced to change
the names of many of the subroutines because most of the ones
with 'kind' in their names are actually returning types.
Modified Paths:
--------------
DART/branches/development/assim_tools/assim_tools_mod.f90
DART/branches/development/assim_tools/assim_tools_mod.html
DART/branches/development/assim_tools/assim_tools_mod.nml
DART/branches/development/cov_cutoff/cov_cutoff_mod.f90
DART/branches/development/location/annulus/location_mod.f90
DART/branches/development/location/channel/location_mod.f90
DART/branches/development/location/column/location_mod.f90
DART/branches/development/location/oned/location_mod.f90
DART/branches/development/location/threed/location_mod.f90
DART/branches/development/location/threed_cartesian/location_mod.f90
DART/branches/development/location/threed_cartesian/xyz_location_mod.f90
DART/branches/development/location/threed_sphere/location_mod.f90
DART/branches/development/location/twod/location_mod.f90
DART/branches/development/location/twod_annulus/location_mod.f90
DART/branches/development/location/twod_sphere/location_mod.f90
DART/branches/development/obs_kind/DEFAULT_obs_kind_mod.F90
-------------- next part --------------
Modified: DART/branches/development/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/branches/development/assim_tools/assim_tools_mod.f90 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/assim_tools/assim_tools_mod.f90 2013-05-30 18:47:39 UTC (rev 6200)
@@ -28,9 +28,9 @@
use obs_def_mod, only : obs_def_type, get_obs_def_location, get_obs_def_time, &
get_obs_def_error_variance, get_obs_kind
-use obs_kind_mod, only : get_num_obs_kinds, get_obs_kind_index
+use obs_kind_mod, only : get_num_obs_kinds, get_obs_kind_index, &
+ get_obs_kind_var_type, assimilate_this_obs_kind
-
use cov_cutoff_mod, only : comp_cov_factor
use reg_factor_mod, only : comp_reg_factor
@@ -41,7 +41,7 @@
use ensemble_manager_mod, only : ensemble_type, get_my_num_vars, get_my_vars, &
compute_copy_mean_var, get_var_owner_index, &
- map_pe_to_task
+ prepare_to_update_copies, map_pe_to_task
use mpi_utilities_mod, only : my_task_id, broadcast_send, broadcast_recv, &
sum_across_tasks
@@ -57,9 +57,6 @@
use assim_model_mod, only : get_state_meta_data, get_close_maxdist_init, &
get_close_obs_init, get_close_obs
-! GSR add routine to check if kind is assimilated ob type
-use obs_kind_mod, only : assimilate_this_obs_kind
-
implicit none
private
@@ -78,7 +75,7 @@
integer :: num_types = 0
real(r8), allocatable :: cutoff_list(:)
logical :: has_special_cutoffs
-character(len = 129) :: errstring
+character(len = 255) :: msgstring, msgstring2, msgstring3
! Need to read in table for off-line based sampling correction and store it
logical :: first_get_correction = .true.
@@ -113,6 +110,7 @@
integer :: adaptive_localization_threshold = -1
real(r8) :: adaptive_cutoff_floor = 0.0_r8
integer :: print_every_nth_obs = 0
+logical :: close_obs_caching = .true.
! since this is in the namelist, it has to have a fixed size.
integer, parameter :: MAX_ITEMS = 300
@@ -148,6 +146,7 @@
integer :: iunit, io, i, j
integer :: num_special_cutoff, type_index
+logical :: namelist_cache_override = .false.
call register_module(source, revision, revdate)
@@ -168,8 +167,8 @@
! FOR NOW, can only do spread restoration with filter option 1 (need to extend this)
if(spread_restoration .and. .not. filter_kind == 1) then
- write(errstring, *) 'cannot combine spread_restoration and filter_kind ', filter_kind
- call error_handler(E_ERR,'assim_tools_init:', errstring, source, revision, revdate)
+ write(msgstring, *) 'cannot combine spread_restoration and filter_kind ', filter_kind
+ call error_handler(E_ERR,'assim_tools_init:', msgstring, source, revision, revdate)
endif
! allocate a list in all cases - even the ones where there is only
@@ -187,12 +186,12 @@
do i = 1, MAX_ITEMS
if(special_localization_obs_types(i) == 'null') exit
if(special_localization_cutoffs(i) == MISSING_R8) then
- write(errstring, *) 'cutoff value', i, ' is uninitialized.'
+ write(msgstring, *) 'cutoff value', i, ' is uninitialized.'
call error_handler(E_ERR,'assim_tools_init:', &
'special cutoff namelist for types and distances do not match', &
source, revision, revdate, &
text2='kind = '//trim(special_localization_obs_types(i)), &
- text3=trim(errstring))
+ text3=trim(msgstring))
endif
j = j + 1
enddo
@@ -203,20 +202,27 @@
do i = 1, num_special_cutoff
type_index = get_obs_kind_index(special_localization_obs_types(i))
if (type_index < 0) then
- write(errstring, *) 'unrecognized TYPE_ in the special localization namelist:'
- call error_handler(E_ERR,'assim_tools_init:', errstring, source, revision, revdate, &
+ write(msgstring, *) 'unrecognized TYPE_ in the special localization namelist:'
+ call error_handler(E_ERR,'assim_tools_init:', msgstring, source, revision, revdate, &
text2=trim(special_localization_obs_types(i)))
endif
cutoff_list(type_index) = special_localization_cutoffs(i)
end do
+! cannot cache previous obs location if different obs types have different
+! localization radii. change it to false, and warn user why.
+if (has_special_cutoffs .and. close_obs_caching) then
+ namelist_cache_override = .true.
+ close_obs_caching = .false.
+endif
+
if (do_output()) then
- write(errstring, '(A,F18.6)') 'The cutoff namelist value is ', cutoff
- call error_handler(E_MSG,'assim_tools_init:', errstring)
- write(errstring, '(A)') 'cutoff is the localization half-width parameter,'
- call error_handler(E_MSG,'assim_tools_init:', errstring)
- write(errstring, '(A,F18.6)') 'so the effective localization radius is ', cutoff*2.0_r8
- call error_handler(E_MSG,'assim_tools_init:', errstring)
+ write(msgstring, '(A,F18.6)') 'The cutoff namelist value is ', cutoff
+ call error_handler(E_MSG,'assim_tools_init:', msgstring)
+ write(msgstring, '(A)') 'cutoff is the localization half-width parameter,'
+ call error_handler(E_MSG,'assim_tools_init:', msgstring)
+ write(msgstring, '(A,F18.6)') 'so the effective localization radius is ', cutoff*2.0_r8
+ call error_handler(E_MSG,'assim_tools_init:', msgstring)
if (has_special_cutoffs) then
call error_handler(E_MSG, '', '')
@@ -224,24 +230,28 @@
call error_handler(E_MSG,'assim_tools_init:','(type name, specified cutoff distance, effective localization radius)')
do i = 1, num_special_cutoff
- write(errstring, '(A32,F18.6,F18.6)') special_localization_obs_types(i), &
+ write(msgstring, '(A32,F18.6,F18.6)') special_localization_obs_types(i), &
special_localization_cutoffs(i), special_localization_cutoffs(i)*2.0_r8
- call error_handler(E_MSG,'assim_tools_init:', errstring)
+ call error_handler(E_MSG,'assim_tools_init:', msgstring)
end do
call error_handler(E_MSG,'assim_tools_init:','all other observation types will use the default cutoff distance')
call error_handler(E_MSG, '', '')
endif
+ if (namelist_cache_override) then
+ call error_handler(E_MSG,'assim_tools_init:','Disabling the close_obs_caching because specialized localization')
+ call error_handler(E_MSG,'assim_tools_init:','distances are enabled. The two cannot be used together.')
+ endif
if(adaptive_localization_threshold > 0) then
- write(errstring, '(A,I10,A)') 'Using adaptive localization, threshold ', &
+ write(msgstring, '(A,I10,A)') 'Using adaptive localization, threshold ', &
adaptive_localization_threshold, ' obs'
- call error_handler(E_MSG,'assim_tools_init:', errstring)
+ call error_handler(E_MSG,'assim_tools_init:', msgstring)
if(adaptive_cutoff_floor > 0.0_r8) then
- write(errstring, '(A,F18.6)') 'Minimum cutoff will not go below ', &
+ write(msgstring, '(A,F18.6)') 'Minimum cutoff will not go below ', &
adaptive_cutoff_floor
call error_handler(E_MSG,'assim_tools_init:', 'Using adaptive localization cutoff floor.', &
- text2=errstring)
+ text2=msgstring)
endif
endif
@@ -292,8 +302,8 @@
real(r8) :: last_close_obs_dist(obs_ens_handle%my_num_vars)
real(r8) :: last_close_state_dist(ens_handle%my_num_vars)
-integer :: my_num_obs, i, j, ispc, owner, owners_index, my_num_state
-integer :: my_obs(obs_ens_handle%my_num_vars), my_state(ens_handle%my_num_vars)
+integer :: my_num_obs, i, j, owner, owners_index, my_num_state
+integer :: my_obs_indx(obs_ens_handle%my_num_vars), my_state_indx(ens_handle%my_num_vars)
integer :: this_obs_key, obs_mean_index, obs_var_index
integer :: grp_beg(num_groups), grp_end(num_groups), grp_size, grp_bot, grp_top, group
integer :: close_obs_ind(obs_ens_handle%my_num_vars)
@@ -302,7 +312,8 @@
integer :: last_close_state_ind(ens_handle%my_num_vars)
integer :: num_close_obs, obs_index, num_close_states, state_index
integer :: total_num_close_obs, last_num_close_obs, last_num_close_states
-integer :: base_obs_kind, my_obs_kind(obs_ens_handle%my_num_vars)
+integer :: base_obs_kind, base_obs_type, my_obs_kind(obs_ens_handle%my_num_vars)
+integer :: my_obs_type(obs_ens_handle%my_num_vars)
integer :: my_state_kind(ens_handle%my_num_vars), nth_obs
integer :: num_close_obs_buffered, num_close_states_buffered
integer :: num_close_obs_calls_made, num_close_states_calls_made
@@ -312,7 +323,7 @@
type(location_type) :: my_obs_loc(obs_ens_handle%my_num_vars)
type(location_type) :: base_obs_loc, last_base_obs_loc, last_base_states_loc
-type(location_type) :: my_state_loc(ens_handle%my_num_vars)
+type(location_type) :: my_state_loc(ens_handle%my_num_vars), dummyloc
type(get_close_type) :: gc_obs, gc_state
type(obs_type) :: observation
type(obs_def_type) :: obs_def
@@ -324,6 +335,10 @@
logical :: local_obs_inflate
logical :: get_close_buffering
+! we are going to read/write the copies array
+call prepare_to_update_copies(ens_handle)
+call prepare_to_update_copies(obs_ens_handle)
+
! Initialize assim_tools_module if needed
if(.not. module_initialized) then
call assim_tools_init
@@ -335,8 +350,13 @@
localization_unit = open_file(localization_diagnostics_file, action = 'append')
endif
-! turn on and off the close buffering
-get_close_buffering = .true.
+! turn on and off the close buffering - you want to turn this off if you
+! have a specialized get_close_obs() routine in your model_mod and you
+! use the type of the observation to determine anything about the distance
+! to the other locations. it will run slower, but will be accurate. if
+! you have not special handling or don't depend on the obs type, leave
+! this with the default .true. value.
+get_close_buffering = close_obs_caching
! For performance, make local copies of these settings which
! are really in the inflate derived type.
@@ -373,7 +393,7 @@
! Get info on my number and indices for obs
my_num_obs = get_my_num_vars(obs_ens_handle)
-call get_my_vars(obs_ens_handle, my_obs)
+call get_my_vars(obs_ens_handle, my_obs_indx)
! Construct an observation temporary
call init_obs(observation, get_num_copies(obs_seq), get_num_qc(obs_seq))
@@ -384,18 +404,23 @@
call get_obs_from_key(obs_seq, this_obs_key, observation)
call get_obs_def(observation, obs_def)
my_obs_loc(i) = get_obs_def_location(obs_def)
- my_obs_kind(i) = get_obs_kind(obs_def)
+ my_obs_type(i) = get_obs_kind(obs_def)
+ if (my_obs_type(i) > 0) then
+ my_obs_kind(i) = get_obs_kind_var_type(my_obs_type(i))
+ else
+ call get_state_meta_data(-1 * my_obs_type(i), dummyloc, my_obs_kind(i)) ! identity obs
+ endif
! Need the time for regression diagnostics potentially; get from first observation
if(i == 1) obs_time = get_obs_def_time(obs_def)
end do Get_Obs_Locations
! Get info on my number and indices for state
my_num_state = get_my_num_vars(ens_handle)
-call get_my_vars(ens_handle, my_state)
+call get_my_vars(ens_handle, my_state_indx)
! Get the location and kind of all my state variables
do i = 1, ens_handle%my_num_vars
- call get_state_meta_data(my_state(i), my_state_loc(i), my_state_kind(i))
+ call get_state_meta_data(my_state_indx(i), my_state_loc(i), my_state_kind(i))
end do
! PAR: MIGHT BE BETTER TO HAVE ONE PE DEDICATED TO COMPUTING
@@ -467,21 +492,26 @@
! to indicate progress is being made and to allow estimates
! of how long the assim will take.
if (nth_obs == 0) then
- write(errstring, '(A,1x,I8,1x,A,I8)') 'Processing observation ', i, &
+ write(msgstring, '(A,1x,I8,1x,A,I8)') 'Processing observation ', i, &
' of ', obs_ens_handle%num_vars
if (print_timestamps == 0) then
- call error_handler(E_MSG,'filter_assim',errstring)
+ call error_handler(E_MSG,'filter_assim',msgstring)
else
- call timestamp(trim(errstring), pos="brief")
+ call timestamp(trim(msgstring), pos="brief")
endif
endif
- ! Every pe has information about obs sequence
+ ! Every pe has information about the global obs sequence
call get_obs_from_key(obs_seq, keys(i), observation)
call get_obs_def(observation, obs_def)
base_obs_loc = get_obs_def_location(obs_def)
obs_err_var = get_obs_def_error_variance(obs_def)
- base_obs_kind = get_obs_kind(obs_def)
+ base_obs_type = get_obs_kind(obs_def)
+ if (base_obs_type > 0) then
+ base_obs_kind = get_obs_kind_var_type(base_obs_type)
+ else
+ call get_state_meta_data(-1 * base_obs_type, dummyloc, base_obs_kind) ! identity obs
+ endif
! Get the value of the observation
call get_obs_values(observation, obs, obs_val_index)
@@ -593,9 +623,13 @@
if (obs_prior_var(group) < 0.0_r8) obs_prior_var(group) = 0.0_r8
end do
- ! Need to get obs density first in case of adaptive localization
+ ! If we are doing adaptive localization then we need to know the number of
+ ! other observations that are within the localization radius. We may need
+ ! to shrink it, and so we need to know this before doing get_close() for the
+ ! state space (even though the state space increments will be computed and
+ ! applied first).
if (.not. get_close_buffering) then
- call get_close_obs(gc_obs, base_obs_loc, base_obs_kind, my_obs_loc, my_obs_kind, &
+ call get_close_obs(gc_obs, base_obs_loc, base_obs_type, my_obs_loc, my_obs_kind, &
num_close_obs, close_obs_ind, close_obs_dist)
else
if (base_obs_loc == last_base_obs_loc) then
@@ -604,7 +638,7 @@
close_obs_dist(:) = last_close_obs_dist(:)
num_close_obs_buffered = num_close_obs_buffered + 1
else
- call get_close_obs(gc_obs, base_obs_loc, base_obs_kind, my_obs_loc, my_obs_kind, &
+ call get_close_obs(gc_obs, base_obs_loc, base_obs_type, my_obs_loc, my_obs_kind, &
num_close_obs, close_obs_ind, close_obs_dist)
last_base_obs_loc = base_obs_loc
last_num_close_obs = num_close_obs
@@ -617,7 +651,7 @@
! set the cutoff default, keep a copy of the original value, and avoid
! looking up the cutoff in a list if the incoming obs is an identity ob
! (and therefore has a negative kind).
- if (base_obs_kind >= 0) then
+ if (base_obs_kind > 0) then
cutoff_orig = cutoff_list(base_obs_kind)
else
cutoff_orig = cutoff
@@ -629,7 +663,7 @@
if(adaptive_localization_threshold > 0) then
! this does a cross-task sum, so all tasks must make this call.
- total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_kind, &
+ total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, &
close_obs_dist, cutoff_rev*2.0_r8)
! Want expected number of close observations to be reduced to some threshold;
@@ -650,7 +684,7 @@
! do it only when diagnostics are requested.
! this does a cross-task sum, so all tasks must make this call.
- rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_kind, &
+ rev_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, &
close_obs_dist, cutoff_rev*2.0_r8)
@@ -679,7 +713,7 @@
! localization radius, set the diag output. this could be large, use carefully.
! this does a cross-task sum, so all tasks must make this call.
- total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_kind, &
+ total_num_close_obs = count_close(num_close_obs, close_obs_ind, my_obs_type, &
close_obs_dist, cutoff_rev*2.0_r8)
if (my_task_id() == 0) then
@@ -696,7 +730,7 @@
! Now everybody updates their close states
! Find state variables on my process that are close to observation being assimilated
if (.not. get_close_buffering) then
- call get_close_obs(gc_state, base_obs_loc, base_obs_kind, my_state_loc, my_state_kind, &
+ call get_close_obs(gc_state, base_obs_loc, base_obs_type, my_state_loc, my_state_kind, &
num_close_states, close_state_ind, close_state_dist)
else
if (base_obs_loc == last_base_states_loc) then
@@ -705,7 +739,7 @@
close_state_dist(:) = last_close_state_dist(:)
num_close_states_buffered = num_close_states_buffered + 1
else
- call get_close_obs(gc_state, base_obs_loc, base_obs_kind, my_state_loc, my_state_kind, &
+ call get_close_obs(gc_state, base_obs_loc, base_obs_type, my_state_loc, my_state_kind, &
num_close_states, close_state_ind, close_state_dist)
last_base_states_loc = base_obs_loc
last_num_close_states = num_close_states
@@ -731,7 +765,7 @@
! Compute the distance and covariance factor
!PAR URGENT: MAKE SURE THIS INDEXING IS CORRECT; SAME FOR OTHER COMP_COV_FACTOR
cov_factor = comp_cov_factor(close_state_dist(j), cutoff_rev, &
- base_obs_loc, base_obs_kind, my_state_loc(state_index), my_state_kind(state_index))
+ base_obs_loc, base_obs_type, my_state_loc(state_index), my_state_kind(state_index))
! If no weight is indicated, no more to do with this state variable
if(cov_factor <= 0.0_r8) cycle STATE_UPDATE
@@ -761,7 +795,7 @@
else
! Pass the time along with the index for possible diagnostic output
! Compute regression factor for this obs-state pair
- reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state(state_index))
+ reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, my_state_indx(state_index))
endif
! The final factor is the minimum of group regression factor and localization cov_factor
@@ -826,11 +860,11 @@
OBS_UPDATE: do j = 1, num_close_obs
obs_index = close_obs_ind(j)
! Only have to update obs that have not yet been used
- if(my_obs(obs_index) > i) then
+ if(my_obs_indx(obs_index) > i) then
! Compute the distance and the covar_factor
cov_factor = comp_cov_factor(close_obs_dist(j), cutoff_rev, &
- base_obs_loc, base_obs_kind, my_obs_loc(obs_index), my_obs_kind(obs_index))
+ base_obs_loc, base_obs_type, my_obs_loc(obs_index), my_obs_kind(obs_index))
if(cov_factor <= 0.0_r8) cycle OBS_UPDATE
! Loop through and update ensemble members in each group
@@ -850,7 +884,7 @@
! Pass the time along with the index for possible diagnostic output
! Compute regression factor for this obs-state pair
! Negative indicates that this is an observation index
- reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, -1*my_obs(obs_index))
+ reg_factor = comp_reg_factor(num_groups, reg_coef, obs_time, i, -1*my_obs_indx(obs_index))
endif
! Final weight is min of group and localization factors
@@ -883,9 +917,9 @@
call get_close_obs_destroy(gc_obs)
! Assure user we have done something
-write(errstring, '(A,I8,A)') &
+write(msgstring, '(A,I8,A)') &
'Processed', obs_ens_handle%num_vars, ' total observations'
-if (print_trace_details >= 0) call error_handler(E_MSG,'filter_assim:',errstring)
+if (print_trace_details >= 0) call error_handler(E_MSG,'filter_assim:',msgstring)
! diagnostics for stats on saving calls by remembering obs at the same location.
! change .true. to .false. in the line below to remove the output completely.
@@ -964,34 +998,59 @@
prior_var = sum((ens - prior_mean)**2) / (ens_size - 1)
endif
-! If both obs_var and prior_var are 0 don't know what to do
-if(obs_var == 0.0_r8 .and. prior_var == 0.0_r8) call error_handler(E_ERR,&
- 'obs_increment', 'Both obs_var and prior_var are zero. This is inconsistent', &
- source, revision, revdate)
+! If both obs_var and prior_var are 0 there is no right thing to do.
+! if obs_var == 0, delta function -- mean becomes obs value with no spread.
+! if prior_var == 0, obs has no effect -- increment is 0
+! if both true, stop with a fatal error.
+if ((obs_var == 0.0_r8) .and. (prior_var == 0.0_r8)) then
-! Call the appropriate filter option to compute increments for ensemble
-if(filter_kind == 1) then
- call obs_increment_eakf(ens, ens_size, prior_mean, prior_var, &
- obs, obs_var, obs_inc, net_a)
-else if(filter_kind == 2) then
- call obs_increment_enkf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
-else if(filter_kind == 3) then
- call obs_increment_kernel(ens, ens_size, obs, obs_var, obs_inc)
-else if(filter_kind == 4) then
- call obs_increment_particle(ens, ens_size, obs, obs_var, obs_inc)
-else if(filter_kind == 5) then
- call obs_increment_ran_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
-else if(filter_kind == 6) then
- call obs_increment_det_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
-else if(filter_kind == 7) then
- call obs_increment_boxcar(ens, ens_size, obs, obs_var, obs_inc, rel_weights)
-else if(filter_kind == 8) then
- call obs_increment_boxcar2(ens, ens_size, prior_mean, prior_var, &
- obs, obs_var, obs_inc)
-else
- call error_handler(E_ERR,'obs_increment', &
- 'Illegal value of filter_kind in assim_tools namelist [1-8 OK]', &
- source, revision, revdate)
+ ! fail if both obs variance and prior spreads are 0.
+ write(msgstring, *) 'Observation value is ', obs, ' ensemble mean value is ', prior_mean
+ write(msgstring2, *) 'The observation has 0.0 error variance, and the ensemble members have 0.0 spread.'
+ write(msgstring3, *) 'These require inconsistent actions and the algorithm cannot continue.'
+ call error_handler(E_ERR, 'obs_increment', msgstring, &
+ source, revision, revdate, text2=msgstring2, text3=msgstring3)
+
+else if (obs_var == 0.0_r8) then
+
+ ! new mean is obs value, so increments are differences between obs
+ ! value and current value. after applying obs, all state will equal obs.
+ obs_inc(:) = obs - ens
+
+else if (prior_var == 0.0_r8) then
+
+ ! if all state values are the same, nothing changes.
+ obs_inc(:) = 0.0_r8
+
+else
+
+ ! Call the appropriate filter option to compute increments for ensemble
+ ! note that at this point we've taken care of the cases where either the
+ ! obs_var or the prior_var is 0, so the individual routines no longer need
+ ! to have code to test for those cases.
+ if(filter_kind == 1) then
+ call obs_increment_eakf(ens, ens_size, prior_mean, prior_var, &
+ obs, obs_var, obs_inc, net_a)
+ else if(filter_kind == 2) then
+ call obs_increment_enkf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
+ else if(filter_kind == 3) then
+ call obs_increment_kernel(ens, ens_size, obs, obs_var, obs_inc)
+ else if(filter_kind == 4) then
+ call obs_increment_particle(ens, ens_size, obs, obs_var, obs_inc)
+ else if(filter_kind == 5) then
+ call obs_increment_ran_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
+ else if(filter_kind == 6) then
+ call obs_increment_det_kf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
+ else if(filter_kind == 7) then
+ call obs_increment_boxcar(ens, ens_size, obs, obs_var, obs_inc, rel_weights)
+ else if(filter_kind == 8) then
+ call obs_increment_boxcar2(ens, ens_size, prior_mean, prior_var, &
+ obs, obs_var, obs_inc)
+ else
+ call error_handler(E_ERR,'obs_increment', &
+ 'Illegal value of filter_kind in assim_tools namelist [1-8 OK]', &
+ source, revision, revdate)
+ endif
endif
! Add in the extra increments if doing observation space covariance inflation
@@ -1002,15 +1061,21 @@
! By doing it here, can take care of both standard non-deterministic updates
! plus non-deterministic obs space covariance inflation. This is expensive, so
! don't use it if it's not needed.
-if(sort_obs_inc) then
- new_val = ens_in + obs_inc
- ! Sorting to make increments as small as possible
- call index_sort(ens_in, ens_index, ens_size)
- call index_sort(new_val, new_index, ens_size)
- do i = 1, ens_size
- obs_inc(ens_index(i)) = new_val(new_index(i)) - ens_in(ens_index(i))
- end do
-end if
+if (sort_obs_inc) then
+ if (filter_kind == 1) then
+ ! sort does nothing and we skip the code
+ ! FIXME: are there other kinds where the sort has no effect?
+ ! add more cases here.
+ else
+ new_val = ens_in + obs_inc
+ ! Sorting to make increments as small as possible
+ call index_sort(ens_in, ens_index, ens_size)
+ call index_sort(new_val, new_index, ens_size)
+ do i = 1, ens_size
+ obs_inc(ens_index(i)) = new_val(new_index(i)) - ens_in(ens_index(i))
+ end do
+ endif
+endif
! Get the net change in spread if obs space inflation was used
if(do_obs_inflate(inflate)) net_a = net_a * sqrt(my_cov_inflate)
@@ -1032,14 +1097,8 @@
real(r8) :: new_mean, var_ratio
! Compute the new mean
-if (obs_var /= 0.0_r8) then
- var_ratio = obs_var / (prior_var + obs_var)
- new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
-! If obs is a delta function, it becomes new value
-else
- var_ratio = 0.0_r8
- new_mean = obs
-endif
+var_ratio = obs_var / (prior_var + obs_var)
+new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
! Compute sd ratio and shift ensemble
a = sqrt(var_ratio)
@@ -1064,15 +1123,9 @@
real(r8) :: temp_mean, temp_var, new_ens(ens_size), new_var
integer :: i
-if (obs_var /= 0.0_r8) then
- var_ratio = obs_var / (prior_var + obs_var)
- new_var = var_ratio * prior_var
- new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
-else
- var_ratio = 0.0_r8
- new_var = var_ratio * prior_var
- new_mean = obs
-endif
+var_ratio = obs_var / (prior_var + obs_var)
+new_var = var_ratio * prior_var
+new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
! This will reproduce exactly for multiple runs with the same task count,
! but WILL NOT reproduce for a different number of MPI tasks.
@@ -1123,21 +1176,9 @@
real(r8) :: new_mean, var_ratio, temp_var, new_ens(ens_size), new_var
integer :: i
-if (obs_var /= 0.0_r8) then
- var_ratio = obs_var / (prior_var + obs_var)
- new_var = var_ratio * prior_var
- new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
-else
- if (prior_var /= 0.0_r8) then
- var_ratio = 0.0_r8
- new_var = var_ratio * prior_var
- new_mean = obs
- else
- call error_handler(E_ERR,'obs_increment_det_kf', &
- 'Both obs_var and prior_var are zero. This is inconsistent', &
- source, revision, revdate)
- endif
-endif
+var_ratio = obs_var / (prior_var + obs_var)
+new_var = var_ratio * prior_var
+new_mean = var_ratio * (prior_mean + prior_var*obs / obs_var)
! Want a symmetric distribution with kurtosis 3 and variance new_var and mean new_mean
if(ens_size /= 20) then
@@ -1221,13 +1262,10 @@
real(r8), intent(in) :: ens(ens_size), obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
-real(r8) :: a, weight(ens_size), rel_weight(ens_size), cum_weight(0:ens_size)
+real(r8) :: weight(ens_size), rel_weight(ens_size), cum_weight(0:ens_size)
real(r8) :: base, frac, new_val(ens_size), weight_sum
integer :: i, j, indx(ens_size)
-! The factor a is not defined for particle filters
-a = -1.0_r8
-
! Begin by computing a weight for each of the prior ensemble members
do i = 1, ens_size
weight(i) = exp(-1.0_r8 * (ens(i) - obs)**2 / (2.0_r8 * obs_var))
@@ -1294,18 +1332,15 @@
real(r8), intent(in) :: ens(ens_size), prior_mean, prior_var, obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
-real(r8) :: a, obs_var_inv, prior_var_inv, new_var, new_mean(ens_size)
+real(r8) :: obs_var_inv, prior_var_inv, new_var, new_mean(ens_size)
! real(r8) :: sx, s_x2
real(r8) :: temp_mean, temp_obs(ens_size)
integer :: i
-! The factor a is not defined for kernel filters
-a = -1.0_r8
-
! Compute mt_rinv_y (obs error normalized by variance)
obs_var_inv = 1.0_r8 / obs_var
+prior_var_inv = 1.0_r8 / prior_var
-prior_var_inv = 1.0_r8 / prior_var
new_var = 1.0_r8 / (prior_var_inv + obs_var_inv)
! If this is first time through, need to initialize the random sequence.
@@ -1414,13 +1449,12 @@
do i = 1, ens_size
unif = random_uniform(inc_ran_seq)
! Figure out which kernel it's in
- do j = 1, ens_size
+ whichk: do j = 1, ens_size
if(unif < cum_frac(j)) then
kernel = j
- goto 10
+ exit whichk
end if
- end do
-10 continue
+ end do whichk
! Next calculate a unit normal in this kernel
norm = random_gaussian(inc_ran_seq, 0.0_r8, sqrt(new_cov))
@@ -1578,9 +1612,9 @@
else if(ens_size < 10000) then
write(correction_file_name, 41) 'final_full.', ens_size
else
- write(errstring,*)'Trying to use ',ens_size,' model states -- too many.'
+ write(msgstring,*)'Trying to use ',ens_size,' model states -- too many.'
call error_handler(E_ERR,'get_correction_from_file','Use less than 10000 ens members.',&
- source,revision,revdate, text2=errstring)
+ source,revision,revdate, text2=msgstring)
11 format(a11, i1)
21 format(a11, i2)
@@ -1590,8 +1624,8 @@
! Make sure that the correction file exists, else an error
if(.not. file_exist(correction_file_name)) then
- write(errstring,*) 'Correction file ', correction_file_name, ' does not exist'
- call error_handler(E_ERR,'get_correction_from_file',errstring,source,revision,revdate)
+ write(msgstring,*) 'Correction file ', correction_file_name, ' does not exist'
+ call error_handler(E_ERR,'get_correction_from_file',msgstring,source,revision,revdate)
endif
! Read in file to get the expected value of the true correlation given the sample
@@ -2066,8 +2100,8 @@
elseif (adj_r2 >= x(j) .and. adj_r2 <= x(j+1)) then
new_ens(i) = adj_r2
else
- errstring = 'Did not get a satisfactory quadratic root'
- call error_handler(E_ERR, 'obs_increment_boxcar2', errstring, &
+ msgstring = 'Did not get a satisfactory quadratic root'
+ call error_handler(E_ERR, 'obs_increment_boxcar2', msgstring, &
source, revision, revdate)
endif
endif
@@ -2436,18 +2470,18 @@
!--------------------------------------------------------------------
-function count_close(num_close, index_list, my_kinds, dist, maxdist)
- integer, intent(in) :: num_close, index_list(:), my_kinds(:)
+function count_close(num_close, index_list, my_types, dist, maxdist)
+ integer, intent(in) :: num_close, index_list(:), my_types(:)
real(r8), intent(in) :: dist(:), maxdist
integer :: count_close
! return the total number of items from the index_list which
-! are kinds which are going to be assimilated, and within distance.
+! are types which are going to be assimilated, and within distance.
! this excludes items on the eval list only, not listed, or
! items too far away. this routine does a global communication
! so if any MPI tasks make this call, all must.
-integer :: k, thiskind, local_count
+integer :: k, thistype, local_count
local_count = 0
do k=1, num_close
@@ -2455,14 +2489,14 @@
! only accept items closer than limit
if (dist(k) > maxdist) cycle
- ! include identity obs, plus kinds on assim list.
+ ! include identity obs, plus types on assim list.
! you have to do the if tests separately because fortran allows
! both parts of an if(a .or. b) test to be eval'd at the same time.
! you'd be using a negative index if it was an identity obs.
- thiskind = my_kinds(index_list(k))
- if (thiskind < 0) then
+ thistype = my_types(index_list(k))
+ if (thistype < 0) then
local_count = local_count + 1
- else if (assimilate_this_obs_kind(thiskind)) then
+ else if (assimilate_this_obs_kind(thistype)) then
local_count = local_count + 1
endif
end do
Modified: DART/branches/development/assim_tools/assim_tools_mod.html
===================================================================
--- DART/branches/development/assim_tools/assim_tools_mod.html 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/assim_tools/assim_tools_mod.html 2013-05-30 18:47:39 UTC (rev 6200)
@@ -295,7 +295,8 @@
spread_restoration, sampling_error_correction, &
adaptive_localization_threshold, adaptive_cutoff_floor, &
output_localization_diagnostics, localization_diagnostics_file, &
-special_localization_obs_types, special_localization_cutoffs
+special_localization_obs_types, special_localization_cutoffs, &
+close_obs_caching
</pre>
</div>
@@ -419,6 +420,15 @@
or to verify progress is being made in cases with suspected problems.
Default: 0 </TD></TR>
+<TR><!--contents--><TD valign=top>close_obs_caching </TD>
+ <!-- type --><TD valign=top>logical </TD>
+ <!--descript--><TD> Leave this set to the default of .TRUE. unless you are using
+ specialized_localization_cutoffs. In that case to get accurate
+ results, set it to .FALSE.. This also needs to be .FALSE. if
+ you have a get_close_obs() routine that uses the types/kinds of
+ the obs to adjust the distances.
+ Default: .TRUE. </TD></TR>
+
</TABLE>
</div>
Modified: DART/branches/development/assim_tools/assim_tools_mod.nml
===================================================================
--- DART/branches/development/assim_tools/assim_tools_mod.nml 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/assim_tools/assim_tools_mod.nml 2013-05-30 18:47:39 UTC (rev 6200)
@@ -11,6 +11,7 @@
print_every_nth_obs = 0,
rectangular_quadrature = .true.,
gaussian_likelihood_tails = .false.,
+ close_obs_caching = .true.
/
# specify these in the same order, the same number of items
Modified: DART/branches/development/cov_cutoff/cov_cutoff_mod.f90
===================================================================
--- DART/branches/development/cov_cutoff/cov_cutoff_mod.f90 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/cov_cutoff/cov_cutoff_mod.f90 2013-05-30 18:47:39 UTC (rev 6200)
@@ -48,7 +48,7 @@
-function comp_cov_factor(z_in, c, obs_loc, obs_kind, target_loc, target_kind, &
+function comp_cov_factor(z_in, c, obs_loc, obs_type, target_loc, target_kind, &
localization_override)
!----------------------------------------------------------------------
! function comp_cov_factor(z_in, c)
@@ -74,7 +74,7 @@
real(r8), intent(in) :: z_in, c
type(location_type), optional, intent(in) :: obs_loc, target_loc
-integer, optional, intent(in) :: obs_kind, target_kind
+integer, optional, intent(in) :: obs_type, target_kind
integer, optional, intent(in) :: localization_override
real(r8) :: comp_cov_factor
Modified: DART/branches/development/location/annulus/location_mod.f90
===================================================================
--- DART/branches/development/location/annulus/location_mod.f90 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/location/annulus/location_mod.f90 2013-05-30 18:47:39 UTC (rev 6200)
@@ -127,14 +127,14 @@
!----------------------------------------------------------------------------
-function get_dist(loc1, loc2, kind1, kind2)
+function get_dist(loc1, loc2, type1, kind2)
! Compute distance between 2 locations. Right now the distance only
! depends on the horizontal. A namelist option might need to be added
! that supports computing a true 3d distance.
type(location_type), intent(in) :: loc1, loc2
-integer, optional, intent(in) :: kind1, kind2
+integer, optional, intent(in) :: type1, kind2
real(r8) :: get_dist
real(r8) :: x1, y1, x2, y2
@@ -729,7 +729,7 @@
!----------------------------------------------------------------------------
-subroutine get_close_obs(gc, base_obs_loc, base_obs_kind, obs, obs_kind, &
+subroutine get_close_obs(gc, base_obs_loc, base_obs_type, obs, obs_kind, &
num_close, close_ind, dist)
! Return how many locations are closer than cutoff distance, along with a
@@ -738,7 +738,7 @@
type(get_close_type), intent(in) :: gc
type(location_type), intent(in) :: base_obs_loc, obs(:)
-integer, intent(in) :: base_obs_kind, obs_kind(:)
+integer, intent(in) :: base_obs_type, obs_kind(:)
integer, intent(out) :: num_close, close_ind(:)
real(r8), optional, intent(out) :: dist(:)
@@ -757,7 +757,7 @@
! Return list of obs that are within maxdist and their distances
num_close = 0
do i = 1, gc%num
- this_dist = get_dist(base_obs_loc, obs(i), base_obs_kind, obs_kind(i))
+ this_dist = get_dist(base_obs_loc, obs(i), base_obs_type, obs_kind(i))
if(this_dist <= gc%maxdist) then
! Add this ob to the list
num_close = num_close + 1
Modified: DART/branches/development/location/channel/location_mod.f90
===================================================================
--- DART/branches/development/location/channel/location_mod.f90 2013-05-30 17:57:31 UTC (rev 6199)
+++ DART/branches/development/location/channel/location_mod.f90 2013-05-30 18:47:39 UTC (rev 6200)
@@ -85,9 +85,9 @@
type box_type
private
- integer, pointer :: obs_box(:) ! (nobs); List of obs indices in boxes
- integer, pointer :: count(:, :, :) ! (nx, ny, nz); # of obs in each box
- integer, pointer :: start(:, :, :) ! (nx, ny, nz); Start of list of obs in this box
+ integer, pointer :: loc_box(:) ! (nloc); List of loc indices in boxes
+ integer, pointer :: count(:, :, :) ! (nx, ny, nz); # of locs in each box
+ integer, pointer :: start(:, :, :) ! (nx, ny, nz); Start of list of locs in this box
real(r8) :: bot_x, top_x ! extents in x, y, z
real(r8) :: bot_y, top_y
real(r8) :: bot_z, top_z
@@ -120,18 +120,13 @@
! If maxdist stays the same, don't need to do box distance calculations
integer :: last_maxdist = -1.0
-!-----------------------------------------------------------------
-! Namelist with default values
-
-! FIXME: give these better defaults? or put in namelist
-! right now they are recomputed based on the setting of nboxes;
-! each is about the cube root of nboxes, but if the bounds are
-! very different distances in each dim, we may want better control
-! over the ratios of each.
integer :: nx = 10
integer :: ny = 10
integer :: nz = 10
+!-----------------------------------------------------------------
+! Namelist with default values
+
logical :: output_box_info = .false.
integer :: print_box_level = 0
! tuning options
@@ -192,7 +187,7 @@
!----------------------------------------------------------------------------
-function get_dist(loc1, loc2, kind1, kind2)
+function get_dist(loc1, loc2, type1, kind2)
! returns the distance between 2 locations
@@ -207,7 +202,7 @@
!
type(location_type), intent(in) :: loc1, loc2
-integer, optional, intent(in) :: kind1, kind2
+integer, optional, intent(in) :: type1, kind2
real(r8) :: get_dist
real(r8) :: x_dif, y_dif, z_dif
@@ -559,7 +554,7 @@
call nc_check(nf90_def_dim(ncid=ncFileID, name='location', len=LocationDims, &
dimid = LocDimID), 'nc_write_location_atts', 'def_dim:location '//trim(fname))
-! Define the observation location variable and attributes
+! Define the location variable and attributes
call nc_check(nf90_def_var(ncid=ncFileID, name='location', xtype=nf90_double, &
dimids=(/ LocDimID, ObsNumDimID /), varid=VarID), &
@@ -656,13 +651,13 @@
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list