[Dart-dev] [6206] DART/branches/development: Enable checks to test for MISSING_R8 values in the state
nancy at ucar.edu
nancy at ucar.edu
Thu May 30 16:31:17 MDT 2013
Revision: 6206
Author: thoar
Date: 2013-05-30 16:31:16 -0600 (Thu, 30 May 2013)
Log Message:
-----------
Enable checks to test for MISSING_R8 values in the state
before computing things like ensemble means/variances, etc.
CLM (and potentially POP) require that MISSING_R8 values CAN
be a part of the DART state vector. The forward observation operators
must still take appropriate action (and both CLM and POP do).
A namelist item for assim_tools_mod has been addded to indicate
whether or not DART should expect MISSING_R8 values in the state vector.
For instance, CLM has multiple snow layers. One ensemble member
may have 3 layers of snow, one may have 4 ... all the state vectors
must be the same shape, so ... a MISSING_R8 is put into all the
unused snow layers. There are CLM variables that indicate how many
snow layers are in use, for example, which is how the forward
observation operators work.
At present, a DEBUG block in update_from_obs_inc() is disabled.
IF enabled, it will check to see if there are any MISSING_R8 values in
the quantities being used. It is expected to be horribly expensive,
as update_from_obs_inc() gets called TWO bazillion times.
Modified Paths:
--------------
DART/branches/development/assim_tools/assim_tools_mod.f90
DART/branches/development/ensemble_manager/ensemble_manager_mod.f90
DART/branches/development/models/clm/work/input.nml
-------------- next part --------------
Modified: DART/branches/development/assim_tools/assim_tools_mod.f90
===================================================================
--- DART/branches/development/assim_tools/assim_tools_mod.f90 2013-05-30 21:09:12 UTC (rev 6205)
+++ DART/branches/development/assim_tools/assim_tools_mod.f90 2013-05-30 22:31:16 UTC (rev 6206)
@@ -124,6 +124,12 @@
logical :: rectangular_quadrature = .true.
logical :: gaussian_likelihood_tails = .false.
+! Some models are allowed to have MISSING_R8 values in the DART state vector.
+! If they are encountered, it is not necessarily a FATAL error.
+! Most of the time, if a MISSING_R8 is encountered, DART should die.
+! CLM and POP (more?) should have allow_missing_in_state = .true.
+logical :: allow_missing_in_state = .false.
+
! Not in the namelist; this var disables the experimental
! linear and spherical case code in the adaptive localization
! sections. to try out the alternatives, set this to .false.
@@ -134,7 +140,8 @@
adaptive_localization_threshold, adaptive_cutoff_floor, &
print_every_nth_obs, rectangular_quadrature, gaussian_likelihood_tails, &
output_localization_diagnostics, localization_diagnostics_file, &
- special_localization_obs_types, special_localization_cutoffs
+ special_localization_obs_types, special_localization_cutoffs, &
+ allow_missing_in_state
!============================================================================
@@ -334,6 +341,7 @@
logical :: local_varying_ss_inflate
logical :: local_obs_inflate
logical :: get_close_buffering
+logical :: missing_in_state
! we are going to read/write the copies array
call prepare_to_update_copies(ens_handle)
@@ -753,6 +761,23 @@
STATE_UPDATE: do j = 1, num_close_states
state_index = close_state_ind(j)
+ ! Some models can take evasive action if one or more of the ensembles have
+ ! a missing value. Generally means 'do nothing' (as opposed to DIE)
+ missing_in_state = any(ens_handle%copies(1:ens_size, state_index) == MISSING_R8)
+
+ if ( missing_in_state ) then
+ if ( allow_missing_in_state ) then
+ cycle STATE_UPDATE
+ else
+ ! FIXME ... at some point ... convey which instances are missing
+ write(msgstring,*)'Encountered a MISSING_R8 in DART at state index ',state_index
+ write(msgstring2,*)'namelist value of allow_missing_in_state (.false.) &
+ &implies a fatal error.'
+ call error_handler(E_ERR, 'filter_assim', msgstring, &
+ source, revision, revdate, text2=msgstring2)
+ endif
+ endif
+
! Get the initial values of inflation for this variable if state varying inflation
if(local_varying_ss_inflate) then
varying_ss_inflate = ens_handle%copies(ENS_INF_COPY, state_index)
@@ -859,9 +884,14 @@
! Now everybody updates their obs priors (only ones after this one)
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_indx(obs_index) > i) then
+ ! If the forward observation operator failed, no need to
+ ! update the unassimilated observations
+ if (any(obs_ens_handle%copies(1:ens_size, obs_index) == MISSING_R8)) cycle OBS_UPDATE
+
! Compute the distance and the covar_factor
cov_factor = comp_cov_factor(close_obs_dist(j), cutoff_rev, &
base_obs_loc, base_obs_type, my_obs_loc(obs_index), my_obs_kind(obs_index))
@@ -1032,7 +1062,7 @@
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)
+ call obs_increment_enkf(ens, ens_size, 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
@@ -1044,8 +1074,7 @@
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)
+ call obs_increment_boxcar2(ens, ens_size, 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]', &
@@ -1321,7 +1350,7 @@
-subroutine obs_increment_enkf(ens, ens_size, prior_mean, prior_var, obs, obs_var, obs_inc)
+subroutine obs_increment_enkf(ens, ens_size, prior_var, obs, obs_var, obs_inc)
!========================================================================
! subroutine obs_increment_enkf(ens, ens_size, obs, obs_var, obs_inc)
!
@@ -1329,7 +1358,7 @@
! ENKF version of obs increment
integer, intent(in) :: ens_size
-real(r8), intent(in) :: ens(ens_size), prior_mean, prior_var, obs, obs_var
+real(r8), intent(in) :: ens(ens_size), prior_var, obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
real(r8) :: obs_var_inv, prior_var_inv, new_var, new_mean(ens_size)
@@ -1488,7 +1517,27 @@
real(r8) :: restoration_inc(ens_size), state_mean, state_var, correl
real(r8) :: factor, exp_true_correl, mean_factor
+logical :: missing_in_state = .false.
+logical :: missing_in_obs = .false.
+logical :: missing_in_incs = .false.
+
+! FIXME if there are some missing values in the state or obs
+! we cannot just include them in the math ... not sure if this
+! routine can be called in these situations ... but ...
+
+if (2 == 1) then ! DEBUG VERBOSE
+ missing_in_state = any(state == MISSING_R8)
+ missing_in_obs = any(obs == MISSING_R8)
+ missing_in_incs = any(obs_inc == MISSING_R8)
+
+ if ( missing_in_state .or. missing_in_obs .or. missing_in_incs ) then
+ write(msgstring,*) 'Should not have missing values at this point'
+ call error_handler(E_ERR,'update_from_obs_inc',msgstring,source,revision,revdate)
+ endif
+endif
+
! For efficiency, just compute regression coefficient here unless correl is needed
+
state_mean = sum(state) / ens_size
obs_state_cov = sum( (state - state_mean) * (obs - obs_prior_mean) ) / (ens_size - 1)
@@ -1865,7 +1914,7 @@
-subroutine obs_increment_boxcar2(ens, ens_size, prior_mean, prior_var, &
+subroutine obs_increment_boxcar2(ens, ens_size, prior_var, &
obs, obs_var, obs_inc)
!------------------------------------------------------------------------
!
@@ -1900,7 +1949,7 @@
! jla at ucar.edu if you are interested in trying it.
integer, intent(in) :: ens_size
-real(r8), intent(in) :: ens(ens_size), prior_mean, prior_var, obs, obs_var
+real(r8), intent(in) :: ens(ens_size), prior_var, obs, obs_var
real(r8), intent(out) :: obs_inc(ens_size)
integer :: i, e_ind(ens_size), lowest_box, j
Modified: DART/branches/development/ensemble_manager/ensemble_manager_mod.f90
===================================================================
--- DART/branches/development/ensemble_manager/ensemble_manager_mod.f90 2013-05-30 21:09:12 UTC (rev 6205)
+++ DART/branches/development/ensemble_manager/ensemble_manager_mod.f90 2013-05-30 22:31:16 UTC (rev 6206)
@@ -356,6 +356,7 @@
! To reproduce for varying pe count, need fixed sequence for each copy
call init_random_seq(random_seq, global_copy_index)
do j = 1, ens_handle%num_vars
+ if (ens_handle%vars(j,i) /= MISSING_R8) &
ens_handle%vars(j, i) = random_gaussian(random_seq, ens_handle%vars(j, i), &
perturbation_amplitude)
end do
@@ -1337,7 +1338,7 @@
if (use_copy2var_send_loop .eqv. .true. ) then
! Switched loop index from receiving_pe to sending_pe
! Aim: to make the commication scale better on Yellowstone, as num_pes >> ens_size
-! For small numbers of tasks (32 or less) the recieving_pe loop may be faster.
+! For small numbers of tasks (32 or less) the receiving_pe loop may be faster.
! Namelist option use_copy2var_send_loop can be used to select which
! communication pattern to use
! Default: use sending_pe loop (use_copy2var_send_loop = .true.)
@@ -1479,9 +1480,13 @@
num_copies = end_copy - start_copy + 1
-do i = 1, ens_handle%my_num_vars
- ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
-end do
+MYLOOP : do i = 1, ens_handle%my_num_vars
+ if (any(ens_handle%copies(start_copy:end_copy, i) == MISSING_R8)) then
+ ens_handle%copies(mean_copy, i) = MISSING_R8
+ else
+ ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
+ endif
+end do MYLOOP
end subroutine compute_copy_mean
@@ -1507,13 +1512,21 @@
num_copies = end_copy - start_copy + 1
-do i = 1, ens_handle%my_num_vars
- ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
- ens_handle%copies(sd_copy, i) = sqrt((sum((ens_handle%copies(start_copy:end_copy, i) - &
- ens_handle%copies(mean_copy, i))**2) / (num_copies - 1)))
-end do
+MYLOOP : do i = 1, ens_handle%my_num_vars
+ if (any(ens_handle%copies(start_copy:end_copy, i) == MISSING_R8)) then
+ ens_handle%copies(mean_copy, i) = MISSING_R8
+ ens_handle%copies( sd_copy, i) = MISSING_R8
+ else
+ ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
+ ens_handle%copies( sd_copy, i) = sqrt((sum((ens_handle%copies(start_copy:end_copy, i) - &
+ ens_handle%copies(mean_copy, i))**2) / (num_copies - 1)))
+ endif
+
+end do MYLOOP
+
end subroutine compute_copy_mean_sd
+
!--------------------------------------------------------------------------------
subroutine compute_copy_mean_var(ens_handle, start_copy, end_copy, mean_copy, var_copy)
@@ -1537,11 +1550,16 @@
num_copies = end_copy - start_copy + 1
-do i = 1, ens_handle%my_num_vars
- ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
- ens_handle%copies(var_copy, i) = (sum((ens_handle%copies(start_copy:end_copy, i) - &
- ens_handle%copies(mean_copy, i))**2) / (num_copies - 1))
-end do
+MYLOOP : do i = 1, ens_handle%my_num_vars
+ if (any(ens_handle%copies(start_copy:end_copy, i) == MISSING_R8)) then
+ ens_handle%copies(mean_copy, i) = MISSING_R8
+ ens_handle%copies( var_copy, i) = MISSING_R8
+ else
+ ens_handle%copies(mean_copy, i) = sum(ens_handle%copies(start_copy:end_copy, i)) / num_copies
+ ens_handle%copies( var_copy, i) = (sum((ens_handle%copies(start_copy:end_copy, i) - &
+ ens_handle%copies(mean_copy, i))**2) / (num_copies - 1))
+ endif
+end do MYLOOP
end subroutine compute_copy_mean_var
Modified: DART/branches/development/models/clm/work/input.nml
===================================================================
--- DART/branches/development/models/clm/work/input.nml 2013-05-30 21:09:12 UTC (rev 6205)
+++ DART/branches/development/models/clm/work/input.nml 2013-05-30 22:31:16 UTC (rev 6206)
@@ -81,6 +81,7 @@
&assim_tools_nml
filter_kind = 1,
cutoff = 1.00,
+ allow_missing_in_state = .true.,
sort_obs_inc = .true.,
spread_restoration = .false.,
sampling_error_correction = .false.,
@@ -171,21 +172,21 @@
'H2OSOI_LIQ', 'KIND_SOIL_MOISTURE',
'H2OSOI_ICE', 'KIND_ICE',
'T_SOISNO', 'KIND_SOIL_TEMPERATURE',
+ 'cpool', 'KIND_CARBON',
+ 'frootc', 'KIND_ROOT_CARBON',
+ 'leafc', 'KIND_LEAF_CARBON',
+ 'leafn', 'KIND_LEAF_NITROGEN',
/
- 'frootc', 'KIND_ROOT_CARBON',
'livecrootc', 'KIND_ROOT_CARBON',
'deadcrootc', 'KIND_ROOT_CARBON',
'livestemc', 'KIND_STEM_CARBON',
'deadstemc', 'KIND_STEM_CARBON',
- 'leafc', 'KIND_LEAF_CARBON',
- 'cpool', 'KIND_CARBON',
'frootn', 'KIND_ROOT_NITROGEN',
'livecrootn', 'KIND_ROOT_NITROGEN',
'deadcrootn', 'KIND_ROOT_NITROGEN',
'livestemn', 'KIND_STEM_NITROGEN',
'deadstemn', 'KIND_STEM_NITROGEN',
- 'leafn', 'KIND_LEAF_NITROGEN',
'litr1c', 'KIND_LEAF_CARBON',
'litr2c', 'KIND_LEAF_CARBON',
'litr3c', 'KIND_LEAF_CARBON',
@@ -227,7 +228,7 @@
&utilities_nml
- TERMLEVEL = 1,
+ TERMLEVEL = 2,
module_details = .false.,
logfilename = 'dart_log.out',
nmlfilename = 'dart_log.nml',
More information about the Dart-dev
mailing list