[Dart-dev] [6170] DART/branches/development: Add routines which verify which is the most recently updated part of the
nancy at ucar.edu
nancy at ucar.edu
Wed May 29 09:54:12 MDT 2013
Revision: 6170
Author: nancy
Date: 2013-05-29 09:54:11 -0600 (Wed, 29 May 2013)
Log Message:
-----------
Add routines which verify which is the most recently updated part of the
ensemble handle before accessing it directly.
Modified Paths:
--------------
DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
DART/branches/development/adaptive_inflate/fill_inflation_restart.f90
DART/branches/development/ensemble_manager/ensemble_manager_mod.f90
DART/branches/development/filter/filter.dopplerfold.f90
DART/branches/development/filter/filter.f90
DART/branches/development/integrate_model/integrate_model.f90
DART/branches/development/obs_model/obs_model_mod.f90
DART/branches/development/perfect_model_obs/perfect_model_obs.f90
DART/branches/development/smoother/smoother_mod.f90
DART/branches/development/utilities/restart_file_tool.f90
-------------- next part --------------
Modified: DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90
===================================================================
--- DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2013-05-28 22:35:54 UTC (rev 6169)
+++ DART/branches/development/adaptive_inflate/adaptive_inflate_mod.f90 2013-05-29 15:54:11 UTC (rev 6170)
@@ -19,7 +19,7 @@
use random_seq_mod, only : random_seq_type, random_gaussian, init_random_seq
use ensemble_manager_mod, only : ensemble_type, all_vars_to_all_copies, all_copies_to_all_vars, &
read_ensemble_restart, write_ensemble_restart, get_copy_owner_index, &
- map_pe_to_task
+ prepare_to_write_to_vars, prepare_to_read_from_vars, prepare_to_update_vars, map_pe_to_task
use mpi_utilities_mod, only : my_task_id, send_to, receive_from
implicit none
@@ -206,11 +206,17 @@
! a transpose.
if (.not. mean_from_restart) then
call get_copy_owner_index(ss_inflate_index, owner, owners_index)
- if (owner == ens_handle%my_pe) ens_handle%vars(:, owners_index) = inf_initial
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_write_to_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = inf_initial
+ endif
endif
if (.not. sd_from_restart) then
call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
- if (owner == ens_handle%my_pe) ens_handle%vars(:, owners_index) = sd_initial
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_write_to_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = sd_initial
+ endif
endif
endif
@@ -225,6 +231,7 @@
! if inflation array is already on PE0, just figure out the
! largest value in the array and we're done.
if (owner == 0) then
+ call prepare_to_read_from_vars(ens_handle)
minmax_mean(1) = minval(ens_handle%vars(:, owners_index))
minmax_mean(2) = maxval(ens_handle%vars(:, owners_index))
else
@@ -234,6 +241,7 @@
if (ens_handle%my_pe == 0) then
call receive_from(map_pe_to_task(ens_handle, owner), minmax_mean)
else if (ens_handle%my_pe == owner) then
+ call prepare_to_read_from_vars(ens_handle)
minmax_mean(1) = minval(ens_handle%vars(:, owners_index))
minmax_mean(2) = maxval(ens_handle%vars(:, owners_index))
call send_to(map_pe_to_task(ens_handle, 0), minmax_mean)
@@ -245,6 +253,7 @@
! if inflation sd array is already on PE0, just figure out the
! largest value in the array and we're done.
if (owner == 0) then
+ call prepare_to_read_from_vars(ens_handle)
minmax_sd(1) = minval(ens_handle%vars(:, owners_index))
minmax_sd(2) = maxval(ens_handle%vars(:, owners_index))
else
@@ -254,6 +263,7 @@
if (ens_handle%my_pe == 0) then
call receive_from(map_pe_to_task(ens_handle, owner), minmax_sd)
else if (ens_handle%my_pe == owner) then
+ call prepare_to_read_from_vars(ens_handle)
minmax_sd(1) = minval(ens_handle%vars(:, owners_index))
minmax_sd(2) = maxval(ens_handle%vars(:, owners_index))
call send_to(map_pe_to_task(ens_handle, 0), minmax_sd)
@@ -269,9 +279,15 @@
! ensure the entire array contains a single constant value to match what the code uses.
if(inf_flavor == 3) then
call get_copy_owner_index(ss_inflate_index, owner, owners_index)
- if (owner == ens_handle%my_pe) ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_update_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ endif
call get_copy_owner_index(ss_inflate_sd_index, owner, owners_index)
- if (owner == ens_handle%my_pe) ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ if (owner == ens_handle%my_pe) then
+ call prepare_to_update_vars(ens_handle)
+ ens_handle%vars(:, owners_index) = ens_handle%vars(1, owners_index)
+ endif
endif
!------ Block for obs. space inflation initialization ------
Modified: DART/branches/development/adaptive_inflate/fill_inflation_restart.f90
===================================================================
--- DART/branches/development/adaptive_inflate/fill_inflation_restart.f90 2013-05-28 22:35:54 UTC (rev 6169)
+++ DART/branches/development/adaptive_inflate/fill_inflation_restart.f90 2013-05-29 15:54:11 UTC (rev 6170)
@@ -20,8 +20,9 @@
use utilities_mod, only : error_handler, E_ERR, E_MSG, &
initialize_utilities, finalize_utilities
-use ensemble_manager_mod, only : ensemble_type, write_ensemble_restart, &
- init_ensemble_manager, end_ensemble_manager
+use ensemble_manager_mod, only : ensemble_type, write_ensemble_restart, &
+ init_ensemble_manager, end_ensemble_manager, &
+ prepare_to_write_to_vars
use assim_model_mod, only : static_init_assim_model
@@ -75,6 +76,7 @@
call error_handler(E_MSG,'',msgstring)
call init_ensemble_manager(ens_handle, 2, model_size)
+call prepare_to_write_to_vars(ens_handle)
ens_handle%vars(:, ss_inflate_index ) = inf_initial
ens_handle%vars(:, ss_inflate_sd_index) = sd_initial
Modified: DART/branches/development/ensemble_manager/ensemble_manager_mod.f90
===================================================================
--- DART/branches/development/ensemble_manager/ensemble_manager_mod.f90 2013-05-28 22:35:54 UTC (rev 6169)
+++ DART/branches/development/ensemble_manager/ensemble_manager_mod.f90 2013-05-29 15:54:11 UTC (rev 6170)
@@ -28,7 +28,8 @@
close_restart, pert_model_state
use time_manager_mod, only : time_type, set_time
use random_seq_mod, only : random_seq_type, init_random_seq, random_gaussian
-use mpi_utilities_mod, only : task_count, my_task_id, send_to, receive_from, task_sync
+use mpi_utilities_mod, only : task_count, my_task_id, send_to, receive_from, &
+ task_sync, broadcast_send, broadcast_recv
use sort_mod, only : index_sort
implicit none
@@ -51,6 +52,9 @@
get_copy, put_copy, all_vars_to_all_copies, &
all_copies_to_all_vars, read_ensemble_restart, write_ensemble_restart, &
compute_copy_mean_var, get_copy_owner_index, set_ensemble_time, &
+ broadcast_copy, prepare_to_write_to_vars, prepare_to_write_to_copies, &
+ prepare_to_read_from_vars, prepare_to_read_from_copies, prepare_to_update_vars, &
+ prepare_to_update_copies, print_ens_handle, &
map_task_to_pe, map_pe_to_task
type ensemble_type
@@ -65,6 +69,8 @@
! Time is only related to var complete
type(time_type), pointer :: time(:)
integer :: distribution_type
+ integer :: valid ! copies modified last, vars modified last, both same
+ integer :: id_num
integer, allocatable :: task_to_pe_list(:), pe_to_task_list(:) ! List of tasks
! Flexible my_pe, layout_type which allows different task layouts for different ensemble handles
integer :: my_pe
@@ -72,15 +78,29 @@
end type ensemble_type
+! track if copies modified last, vars modified last, both are in sync
+! (and therefore both valid to be used r/o), or unknown.
+integer, parameter :: VALID_UNKNOWN = -1
+integer, parameter :: VALID_BOTH = 0 ! vars & copies have same data
+integer, parameter :: VALID_VARS = 1 ! vars(:,:) modified last
+integer, parameter :: VALID_COPIES = 2 ! copies(:,:) modified last
+
+! unique counter per ensemble handle
+integer :: global_counter = 1
+
! Logical flag for initialization of module
logical :: module_initialized = .false.
! Module storage for writing error messages
-character(len = 129) :: errstring
+character(len = 255) :: msgstring
! Module storage for pe information for this process avoids recomputation
integer :: num_pes
+! Control order of communication loops in the transpose routines
+logical :: use_copy2var_send_loop = .true.
+logical :: use_var2copy_rec_loop = .true.
+
!-----------------------------------------------------------------
!
! namelist with default values
@@ -93,9 +113,11 @@
logical :: single_restart_file_out = .true.
! Size of perturbations for creating ensembles when model won't do it
real(r8) :: perturbation_amplitude = 0.2_r8
-! Options to change order of communication loops in the transpose routines
-logical :: use_copy2var_send_loop = .true.
-logical :: use_var2copy_rec_loop = .true.
+! Complain if unneeded transposes are done
+logical :: flag_unneeded_transposes = .false.
+! Communication configuration:
+! 1 = usual default, 2 - 4 are valid and depend on the machine, ensemble count, and task count
+integer :: communication_configuration = 1
! task layout options:
integer :: layout = 1 ! default to my_pe = my_task_id(). Layout2 assumes that the user knows the correct tasks_per_node
integer :: tasks_per_node = 1 ! default to 1 if the user does not specify a number of tasks per node.
@@ -104,10 +126,9 @@
namelist / ensemble_manager_nml / single_restart_file_in, &
single_restart_file_out, &
perturbation_amplitude, &
- use_copy2var_send_loop, &
- use_var2copy_rec_loop, &
+ communication_configuration, &
layout, tasks_per_node, &
- debug
+ debug, flag_unneeded_transposes
!-----------------------------------------------------------------
@@ -164,7 +185,8 @@
! Check for error: only layout_types 1 and 2 are implemented
if (ens_handle%layout_type /= 1 .and. ens_handle%layout_type /=2 ) then
- call error_handler(E_ERR, 'init_ensemble_manager', 'only layout values 1 (standard), 2(round-robin) allowed ', source, revision, revdate)
+ call error_handler(E_ERR, 'init_ensemble_manager', 'only layout values 1 (standard), 2 (round-robin) allowed ', &
+ source, revision, revdate)
endif
allocate(ens_handle%task_to_pe_list(num_pes))
@@ -177,6 +199,37 @@
ens_handle%num_copies = num_copies
ens_handle%num_vars = num_vars
+! For debugging, error checking
+ens_handle%id_num = global_counter
+global_counter = global_counter + 1
+
+! This controls which way the two transpose routines order their
+! loops: single sender with multiple receivers, or single receiver
+! with multiple senders. For small problems it doesn't matter;
+! for very large ensemble sizes and very large MPI task counts you
+! will want to profile each combination and pick the fastest.
+if (communication_configuration == 1) then
+ use_copy2var_send_loop = .true.
+ use_var2copy_rec_loop = .true.
+else if (communication_configuration == 2) then
+ use_copy2var_send_loop = .false.
+ use_var2copy_rec_loop = .true.
+else if (communication_configuration == 3) then
+ use_copy2var_send_loop = .true.
+ use_var2copy_rec_loop = .false.
+else if (communication_configuration == 4) then
+ use_copy2var_send_loop = .false.
+ use_var2copy_rec_loop = .false.
+else
+ write(msgstring, *) 'communication_configuration is', communication_configuration
+ call error_handler(E_ERR, 'init_ensemble_manager', &
+ 'communication_configuration must be between 1 and 4', &
+ source, revision, revdate, text2=msgstring)
+endif
+
+! initially no data
+ens_handle%valid = VALID_UNKNOWN
+
if(debug .and. my_task_id()==0) then
print*, 'pe_to_task_list', ens_handle%pe_to_task_list
print*, 'task_to_pe_list', ens_handle%task_to_pe_list
@@ -185,6 +238,8 @@
! Figure out how the ensemble copies are partitioned
call set_up_ens_distribution(ens_handle)
+if (debug) call print_ens_handle(ens_handle)
+
end subroutine init_ensemble_manager
!-------------------------------------------------------------------------------
@@ -217,11 +272,13 @@
! Does not make sense to have start_from_restart and single_restart_file_in BOTH false
if(.not. start_from_restart .and. .not. single_restart_file_in) then
- write(errstring, *) 'start_from_restart in filter_nml and single_restart_file_in in &
+ write(msgstring, *) 'start_from_restart in filter_nml and single_restart_file_in in &
&ensemble_manager_nml cannot both be false'
- call error_handler(E_ERR,'read_ensemble_restart', errstring, source, revision, revdate)
+ call error_handler(E_ERR,'read_ensemble_restart', msgstring, source, revision, revdate)
endif
+! the code reads into the vars array
+ens_handle%valid = VALID_VARS
! Some compilers (absoft, but others also) are particularly unhappy about
! both checking present(N) _and_ evaluating N inside a single if() test.
@@ -334,6 +391,12 @@
single_file_forced = .FALSE.
endif
+! Error checking
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'write_ensemble_restart', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+
! For single file, need to send restarts to pe0 and it writes them out.
!-------------- Block for single_restart file -------------
! Need to force single restart file for inflation files
@@ -416,16 +479,22 @@
integer :: owner, owners_index
+! Error checking
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'get_copy', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+
! Verify that requested copy exists
if(copy < 1 .or. copy > ens_handle%num_copies) then
- write(errstring, *) 'Requested copy: ', copy, ' is > maximum copy: ', ens_handle%num_copies
- call error_handler(E_ERR,'get_copy', errstring, source, revision, revdate)
+ write(msgstring, *) 'Requested copy: ', copy, ' is > maximum copy: ', ens_handle%num_copies
+ call error_handler(E_ERR,'get_copy', msgstring, source, revision, revdate)
endif
! Make sure that vars has enough space to handle the answer
if(size(vars) < ens_handle%num_vars) then
- write(errstring, *) 'Size of vars: ', size(vars), ' Must be at least ', ens_handle%num_vars
- call error_handler(E_ERR,'get_copy', errstring, source, revision, revdate)
+ write(msgstring, *) 'Size of vars: ', size(vars), ' Must be at least ', ens_handle%num_vars
+ call error_handler(E_ERR,'get_copy', msgstring, source, revision, revdate)
endif
! Figure out which PE stores this copy and what its local storage index is
@@ -479,15 +548,21 @@
integer :: owner, owners_index
+! Error checking
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'put_copy', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+
if(copy < 1 .or. copy > ens_handle%num_copies) then
- write(errstring, *) 'Requested copy: ', copy, ' is > maximum copy: ', ens_handle%num_copies
- call error_handler(E_ERR,'put_copy', errstring, source, revision, revdate)
+ write(msgstring, *) 'Requested copy: ', copy, ' is > maximum copy: ', ens_handle%num_copies
+ call error_handler(E_ERR,'put_copy', msgstring, source, revision, revdate)
endif
! Make sure that vars has enough space to handle the answer
if(size(vars) < ens_handle%num_vars) then
- write(errstring, *) 'Size of vars: ', size(vars), ' Must be at least ', ens_handle%num_vars
- call error_handler(E_ERR,'put_copy', errstring, source, revision, revdate)
+ write(msgstring, *) 'Size of vars: ', size(vars), ' Must be at least ', ens_handle%num_vars
+ call error_handler(E_ERR,'put_copy', msgstring, source, revision, revdate)
endif
! What PE stores this copy and what is its local storage index
@@ -522,6 +597,139 @@
!-----------------------------------------------------------------
+subroutine broadcast_copy(ens_handle, copy, arraydata)
+
+! find which PE has the global copy number and have it broadcast
+! that copy to all the other PEs. arraydata is an output on
+! all PEs, even on the PE which is the owner it is separate
+! storage from the vars array in the ensemble handle.
+
+type(ensemble_type), intent(in) :: ens_handle
+integer, intent(in) :: copy
+real(r8), intent(out) :: arraydata(:)
+
+integer :: owner, owners_index
+
+! Error checking
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'broadcast_copy', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+
+if(copy < 1 .or. copy > ens_handle%num_copies) then
+ write(msgstring, *) 'Requested copy: ', copy, ' is > maximum copy: ', ens_handle%num_copies
+ call error_handler(E_ERR,'broadcast_copy', msgstring, source, revision, revdate)
+endif
+
+! Make sure that arraydata has enough space to handle the answer
+if(size(arraydata) < ens_handle%num_vars) then
+ write(msgstring, *) 'Size of arraydata: ', size(arraydata), ' must be at least ', ens_handle%num_vars
+ call error_handler(E_ERR,'broadcast_copy', msgstring, source, revision, revdate)
+endif
+
+! What PE stores this copy and what is its local storage index
+call get_copy_owner_index(copy, owner, owners_index)
+
+! First block of code that must be done by PE that is to send the copy
+if(ens_handle%my_pe == owner) then
+ arraydata = ens_handle%vars(:, owners_index)
+ call broadcast_send(map_pe_to_task(ens_handle, owner), arraydata)
+else
+ call broadcast_recv(map_pe_to_task(ens_handle, owner), arraydata)
+endif
+
+end subroutine broadcast_copy
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_write_to_vars(ens_handle)
+
+! Warn ens manager that we're going to directly update the %vars array
+
+type(ensemble_type), intent(inout) :: ens_handle
+
+ens_handle%valid = VALID_VARS
+
+end subroutine prepare_to_write_to_vars
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_write_to_copies(ens_handle)
+
+! Warn ens manager that we're going to directly update the %copies array
+
+type(ensemble_type), intent(inout) :: ens_handle
+
+ens_handle%valid = VALID_COPIES
+
+end subroutine prepare_to_write_to_copies
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_read_from_vars(ens_handle)
+
+! Check to be sure that the vars array is current
+
+type(ensemble_type), intent(in) :: ens_handle
+
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'prepare_to_read_from_vars', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+
+end subroutine prepare_to_read_from_vars
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_read_from_copies(ens_handle)
+
+! Check to be sure that the copies array is current
+
+type(ensemble_type), intent(in) :: ens_handle
+
+if (ens_handle%valid /= VALID_COPIES .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'prepare_to_read_from_copies', &
+ 'last access not copy-complete', source, revision, revdate)
+endif
+
+end subroutine prepare_to_read_from_copies
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_update_vars(ens_handle)
+
+! We need read/write access, so it has to start valid for vars or both,
+! and then is going to be vars only going out.
+
+type(ensemble_type), intent(inout) :: ens_handle
+
+if (ens_handle%valid /= VALID_VARS .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'prepare_to_update_vars', &
+ 'last access not var-complete', source, revision, revdate)
+endif
+ens_handle%valid = VALID_VARS
+
+end subroutine prepare_to_update_vars
+
+!-----------------------------------------------------------------
+
+subroutine prepare_to_update_copies(ens_handle)
+
+! We need read/write access, so it has to start valid for copies or both,
+! and then is going to be copies only going out.
+
+type(ensemble_type), intent(inout) :: ens_handle
+
+if (ens_handle%valid /= VALID_COPIES .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'prepare_to_update_copies', &
+ 'last access not copy-complete', source, revision, revdate)
+endif
+ens_handle%valid = VALID_COPIES
+
+end subroutine prepare_to_update_copies
+
+!-----------------------------------------------------------------
+
subroutine set_ensemble_time(ens_handle, indx, mtime)
! Sets the time of an ensemble member indexed by local storage on this pe.
@@ -531,8 +739,8 @@
type(time_type), intent(in) :: mtime
if(indx < 1 .or. indx > ens_handle%my_num_copies) then
- write(errstring, *) 'indx: ', indx, ' cannot exceed ', ens_handle%my_num_copies
- call error_handler(E_ERR,'get_ensemble_time', errstring, source, revision, revdate)
+ write(msgstring, *) 'indx: ', indx, ' cannot exceed ', ens_handle%my_num_copies
+ call error_handler(E_ERR,'get_ensemble_time', msgstring, source, revision, revdate)
endif
ens_handle%time(indx) = mtime
@@ -550,8 +758,8 @@
type(time_type), intent(out) :: mtime
if(indx < 1 .or. indx > ens_handle%my_num_copies) then
- write(errstring, *) 'indx: ', indx, ' cannot exceed ', ens_handle%my_num_copies
- call error_handler(E_ERR,'get_ensemble_time', errstring, source, revision, revdate)
+ write(msgstring, *) 'indx: ', indx, ' cannot exceed ', ens_handle%my_num_copies
+ call error_handler(E_ERR,'get_ensemble_time', msgstring, source, revision, revdate)
endif
mtime = ens_handle%time(indx)
@@ -581,31 +789,39 @@
logical, intent(in) :: duplicate_time
! Name was changed from copy_ens to avoid confusion in naming.
-! Should only be used if the ens%vars storage is current (each pe has subset of
+! Should only be used if the ens1%vars storage is current (each pe has subset of
! copies of all vars).
! If duplicate_time is true, also copies the time information from ens1 to ens2.
! If duplicate_time is false, the times in ens2 are left unchanged.
+! Error checking
+if (ens1%valid /= VALID_VARS .and. ens1%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'duplicate_ens', &
+ 'last access not var-complete for source ensemble', source, revision, revdate)
+endif
+
! Check to make sure that the ensembles are compatible
if(ens1%num_copies /= ens2%num_copies) then
- write(errstring, *) 'num_copies ', ens1%num_copies, ' and ', ens2%num_copies, &
+ write(msgstring, *) 'num_copies ', ens1%num_copies, ' and ', ens2%num_copies, &
'must be equal'
- call error_handler(E_ERR,'duplicate_ens', errstring, source, revision, revdate)
+ call error_handler(E_ERR,'duplicate_ens', msgstring, source, revision, revdate)
endif
if(ens1%num_vars /= ens2%num_vars) then
- write(errstring, *) 'num_vars ', ens1%num_vars, ' and ', ens2%num_vars, &
+ write(msgstring, *) 'num_vars ', ens1%num_vars, ' and ', ens2%num_vars, &
'must be equal'
- call error_handler(E_ERR,'duplicate_ens', errstring, source, revision, revdate)
+ call error_handler(E_ERR,'duplicate_ens', msgstring, source, revision, revdate)
endif
if(ens1%distribution_type /= ens2%distribution_type) then
- write(errstring, *) 'distribution_type ', ens1%distribution_type, ' and ', &
+ write(msgstring, *) 'distribution_type ', ens1%distribution_type, ' and ', &
ens2%distribution_type, 'must be equal'
- call error_handler(E_ERR,'duplicate_ens', errstring, source, revision, revdate)
+ call error_handler(E_ERR,'duplicate_ens', msgstring, source, revision, revdate)
endif
! Duplicate each copy that is stored locally on this process.
ens2%vars = ens1%vars
+ens2%valid = VALID_VARS
+
! Duplicate time if requested
if(duplicate_time) ens2%time = ens1%time
@@ -637,9 +853,9 @@
integer, intent(out) :: copies(:)
if(size(copies) < ens_handle%my_num_copies) then
- write(errstring, *) 'Array copies only has size ', size(copies), &
+ write(msgstring, *) 'Array copies only has size ', size(copies), &
' but must be at least ', ens_handle%my_num_copies
- call error_handler(E_ERR, 'get_my_copies', errstring, source, revision, revdate)
+ call error_handler(E_ERR, 'get_my_copies', msgstring, source, revision, revdate)
endif
copies(1:ens_handle%my_num_copies) = ens_handle%my_copies(1:ens_handle%my_num_copies)
@@ -672,9 +888,9 @@
integer, intent(out) :: vars(:)
if(size(vars) < ens_handle%my_num_vars) then
- write(errstring, *) 'Array vars only has size ', size(vars), &
+ write(msgstring, *) 'Array vars only has size ', size(vars), &
' but must be at least ', ens_handle%my_num_vars
- call error_handler(E_ERR,'get_my_vars', errstring, source, revision, revdate)
+ call error_handler(E_ERR,'get_my_vars', msgstring, source, revision, revdate)
endif
vars(1:ens_handle%my_num_vars) = ens_handle%my_vars(1:ens_handle%my_num_vars)
@@ -905,6 +1121,24 @@
call timestamp_message('vars_to_copies start: '//label, alltasks=.true.)
endif
+! Error checking, and early return if possible.
+if (ens_handle%valid == VALID_BOTH) then
+ if (flag_unneeded_transposes) then
+ write(msgstring, *) 'ens_handle ', ens_handle%id_num
+ call error_handler(E_MSG, 'all_vars_to_all_copies', &
+ 'vars & copies both valid, transpose not needed', source, revision, revdate, &
+ text2=msgstring)
+ endif
+ return
+else if (ens_handle%valid /= VALID_VARS) then
+ write(msgstring, *) 'ens_handle ', ens_handle%id_num
+ call error_handler(E_ERR, 'all_vars_to_all_copies', &
+ 'last access not var-complete', source, revision, revdate, &
+ text2=msgstring)
+endif
+
+ens_handle%valid = VALID_BOTH
+
! Accelerated version for single process
if(num_pes == 1) then
ens_handle%copies = transpose(ens_handle%vars)
@@ -1059,6 +1293,24 @@
call timestamp_message('copies_to_vars start: '//label, alltasks=.true.)
endif
+! Error checking, and early return if possible
+if (ens_handle%valid == VALID_BOTH) then
+ if (flag_unneeded_transposes) then
+ write(msgstring, *) 'ens_handle ', ens_handle%id_num
+ call error_handler(E_MSG, 'all_copies_to_all_vars', &
+ 'vars & copies both valid, transpose not needed', source, revision, revdate, &
+ text2=msgstring)
+ endif
+ return
+else if (ens_handle%valid /= VALID_COPIES) then
+ write(msgstring, *) 'ens_handle ', ens_handle%id_num
+ call error_handler(E_ERR, 'all_copies_to_all_vars', &
+ 'last access not copy-complete', source, revision, revdate, &
+ text2=msgstring)
+endif
+
+ens_handle%valid = VALID_BOTH
+
! Accelerated version for single process
if(num_pes == 1) then
ens_handle%vars = transpose(ens_handle%copies)
@@ -1219,6 +1471,12 @@
! Should check to make sure that start, end and mean are all legal
+! Error checking
+if (ens_handle%valid /= VALID_COPIES .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'compute_copy_mean', &
+ 'last access not copy-complete', source, revision, revdate)
+endif
+
num_copies = end_copy - start_copy + 1
do i = 1, ens_handle%my_num_vars
@@ -1241,6 +1499,12 @@
! Should check to make sure that start, end, mean and sd are all legal copies
+! Error checking
+if (ens_handle%valid /= VALID_COPIES .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'compute_copy_mean_sd', &
+ 'last access not copy-complete', source, revision, revdate)
+endif
+
num_copies = end_copy - start_copy + 1
do i = 1, ens_handle%my_num_vars
@@ -1265,6 +1529,12 @@
! Should check to make sure that start, end, mean and var are all legal copies
+! Error checking
+if (ens_handle%valid /= VALID_COPIES .and. ens_handle%valid /= VALID_BOTH) then
+ call error_handler(E_ERR, 'compute_copy_mean_var', &
+ 'last access not copy-complete', source, revision, revdate)
+endif
+
num_copies = end_copy - start_copy + 1
do i = 1, ens_handle%my_num_vars
@@ -1313,6 +1583,58 @@
!--------------------------------------------------------------------------------
+subroutine print_ens_handle(ens_handle, force, label)
+ type(ensemble_type), intent(in) :: ens_handle
+ logical, optional, intent(in) :: force
+ character(len=*), optional, intent(in) :: label
+
+logical :: print_anyway
+logical :: has_label
+
+print_anyway = .false.
+if (present(force)) then
+ print_anyway = force
+endif
+
+has_label = .false.
+if (present(label)) then
+ has_label = .true.
+endif
+
+! print out contents of an ensemble handle derived type
+if (.not. debug .and. .not. print_anyway) return
+
+if (has_label) then
+ call error_handler(E_MSG, 'ensemble handle: ', label, source, revision, revdate)
+endif
+write(msgstring, *) 'handle num: ', ens_handle%id_num
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'number of copies: ', ens_handle%num_copies
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'number of vars : ', ens_handle%num_vars
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'number of my_copies: ', ens_handle%my_num_copies
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'number of my_vars : ', ens_handle%my_num_vars
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'valid : ', ens_handle%valid
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'distribution_type : ', ens_handle%distribution_type
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+write(msgstring, *) 'my_pe number : ', ens_handle%my_pe
+call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+if (allocated(ens_handle%pe_to_task_list)) then
+ write(msgstring, *) 'task_to_pe_list : ', ens_handle%task_to_pe_list
+ call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+ write(msgstring, *) 'pe_to_task_list : ', ens_handle%pe_to_task_list
+ call error_handler(E_MSG, 'ensemble handle: ', msgstring, source, revision, revdate)
+endif
+
+
+end subroutine print_ens_handle
+
+!--------------------------------------------------------------------------------
+
subroutine assign_tasks_to_pes(ens_handle, nEns_members, layout_type)
! Calulate the task layout based on the tasks per node and the total number of tasks.
! Allows the user to spread out the ensemble members as much as possible to balance
Modified: DART/branches/development/filter/filter.dopplerfold.f90
===================================================================
--- DART/branches/development/filter/filter.dopplerfold.f90 2013-05-28 22:35:54 UTC (rev 6169)
+++ DART/branches/development/filter/filter.dopplerfold.f90 2013-05-29 15:54:11 UTC (rev 6170)
@@ -21,8 +21,9 @@
set_qc_meta_data, get_expected_obs, get_first_obs, &
get_obs_time_range, delete_obs_from_seq, delete_seq_head, &
delete_seq_tail, replace_obs_values, replace_qc, &
- destroy_obs_sequence, get_qc_meta_data
-use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, get_obs_def_time
+ destroy_obs_sequence, get_qc_meta_data, add_qc
+use obs_def_mod, only : obs_def_type, get_obs_def_error_variance, get_obs_def_time, &
+ get_obs_kind
use time_manager_mod, only : time_type, get_time, set_time, operator(/=), operator(>), &
operator(-), print_time
use utilities_mod, only : register_module, error_handler, E_ERR, E_MSG, E_DBG, &
@@ -40,8 +41,10 @@
read_ensemble_restart, write_ensemble_restart, &
compute_copy_mean, compute_copy_mean_sd, &
compute_copy_mean_var, duplicate_ens, get_copy_owner_index, &
- get_ensemble_time, set_ensemble_time, &
- map_task_to_pe, map_pe_to_task
+ get_ensemble_time, set_ensemble_time, broadcast_copy, &
+ prepare_to_read_from_vars, prepare_to_write_to_vars, prepare_to_read_from_copies, &
+ prepare_to_write_to_copies, get_ensemble_time, set_ensemble_time, &
+ map_task_to_pe, map_pe_to_task, prepare_to_update_copies
use adaptive_inflate_mod, only : adaptive_inflate_end, do_varying_ss_inflate, &
do_single_ss_inflate, inflate_ens, adaptive_inflate_init, &
do_obs_inflate, adaptive_inflate_type, &
@@ -276,6 +279,7 @@
call trace_message('After setting up space for observations')
call trace_message('Before setting up space for ensembles')
+
! Allocate model size storage and ens_size storage for metadata for outputting ensembles
model_size = get_model_size()
@@ -509,6 +513,7 @@
if(do_single_ss_inflate(prior_inflate) .or. do_varying_ss_inflate(prior_inflate)) then
call trace_message('Before prior inflation damping and prep')
if (inf_damping(1) /= 1.0_r8) then
+ call prepare_to_update_copies(ens_handle)
ens_handle%copies(PRIOR_INF_COPY, :) = 1.0_r8 + &
inf_damping(1) * (ens_handle%copies(PRIOR_INF_COPY, :) - 1.0_r8)
endif
@@ -531,26 +536,27 @@
! and obs_values. ens_size is the number of regular ensemble members,
! not the number of copies
call get_obs_ens(ens_handle, obs_ens_handle, forward_op_ens_handle, &
- seq, keys, obs_val_index, num_obs_in_set, &
- OBS_ERR_VAR_COPY, OBS_VAL_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY)
+ seq, keys, obs_val_index, input_qc_index, num_obs_in_set, &
+ OBS_ERR_VAR_COPY, OBS_VAL_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
+ isprior=.true.)
! Although they are integer, keys are one 'copy' of obs ensemble
! (the last one?)
call put_copy(map_task_to_pe(obs_ens_handle, 0), obs_ens_handle, OBS_KEY_COPY, keys * 1.0_r8)
- ! Ship the ensemble mean to the model; some models need this for
- ! computing distances. find out who stores the ensemble mean copy.
+
+ ! While we're here, make sure the timestamp on the actual ensemble copy
+ ! for the mean has the current time. If the user requests it be written
+ ! out, it needs a valid timestamp.
call get_copy_owner_index(ENS_MEAN_COPY, mean_owner, mean_owners_index)
- ! Broadcast it to everybody else
if(ens_handle%my_pe == mean_owner) then
! Make sure the timestamp for the mean is the current time.
call set_ensemble_time(ens_handle, mean_owners_index, curr_ens_time)
- ens_mean = ens_handle%vars(:, mean_owners_index)
- call broadcast_send(my_task_id(), ens_mean)
- else
- call broadcast_recv(map_pe_to_task(ens_handle, mean_owner), ens_mean)
endif
- ! Now send the mean to the model in case it's needed
+ ! Make sure all tasks have a copy of the ensemble mean.
+ call broadcast_copy(ens_handle, ENS_MEAN_COPY, ens_mean)
+
+ ! Now send the mean to the model_mod in case it's needed
call ens_mean_for_model(ens_mean)
call timestamp_message('After computing prior observation values')
@@ -581,8 +587,9 @@
OBS_VAL_COPY, OBS_ERR_VAR_COPY, DART_qc_index)
call trace_message('After observation space diagnostics')
- ! Need obs to be copy complete for assimilation
- call all_vars_to_all_copies(obs_ens_handle)
+ ! FIXME: i believe both copies and vars are equal at the end
+ ! of the obs_space diags, so we can skip this.
+ !call all_vars_to_all_copies(obs_ens_handle)
write(msgstring, '(A,I8,A)') 'Ready to assimilate up to', size(keys), ' observations'
call trace_message(msgstring, 'filter:', -1)
@@ -626,6 +633,7 @@
if(do_single_ss_inflate(post_inflate) .or. do_varying_ss_inflate(post_inflate)) then
call trace_message('Before posterior inflation damping and prep')
if (inf_damping(2) /= 1.0_r8) then
+ call prepare_to_update_copies(ens_handle)
ens_handle%copies(POST_INF_COPY, :) = 1.0_r8 + &
inf_damping(2) * (ens_handle%copies(POST_INF_COPY, :) - 1.0_r8)
endif
@@ -650,8 +658,9 @@
! and obs_values. ens_size is the number of regular ensemble members,
! not the number of copies
call get_obs_ens(ens_handle, obs_ens_handle, forward_op_ens_handle, &
- seq, keys, obs_val_index, num_obs_in_set, &
- OBS_ERR_VAR_COPY, OBS_VAL_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY)
+ seq, keys, obs_val_index, input_qc_index, num_obs_in_set, &
+ OBS_ERR_VAR_COPY, OBS_VAL_COPY, OBS_KEY_COPY, OBS_GLOBAL_QC_COPY, &
+ isprior=.false.)
call timestamp_message('After computing posterior observation values')
call trace_message('After computing posterior observation values')
@@ -707,15 +716,7 @@
call timestamp_message('Before computing posterior state space inflation')
! Ship the ensemble mean to the model; some models need this for computing distances
- ! Who stores the ensemble mean copy
- call get_copy_owner_index(ENS_MEAN_COPY, mean_owner, mean_owners_index)
- ! Broadcast it to everybody else
- if(ens_handle%my_pe == mean_owner) then
- ens_mean = ens_handle%vars(:, mean_owners_index)
- call broadcast_send(map_pe_to_task(ens_handle, mean_owner), ens_mean)
- else
- call broadcast_recv(map_pe_to_task(ens_handle, mean_owner), ens_mean)
- endif
+ call broadcast_copy(ens_handle, ENS_MEAN_COPY, ens_mean)
! Now send the mean to the model in case it's needed
call ens_mean_for_model(ens_mean)
@@ -969,8 +970,8 @@
integer, intent(out) :: in_obs_copy, obs_val_index
integer, intent(out) :: input_qc_index, DART_qc_index
-character(len = metadatalength) :: qc_meta_data = 'DART quality control'
character(len = metadatalength) :: no_qc_meta_data = 'No incoming data QC'
+character(len = metadatalength) :: dqc_meta_data = 'DART quality control'
character(len = 129) :: obs_seq_read_format
integer :: obs_seq_file_id, num_obs_copies
integer :: tnum_copies, tnum_qc, tnum_obs, tmax_num_obs, qc_num_inc, num_qc
@@ -981,32 +982,62 @@
! Prior and posterior values for all selected fields (so times 2)
num_obs_copies = 2 * num_output_obs_members + 4
-! Input file can have one qc field or none, not prepared to have more
-! The one that exists would be the NCEP and perfect_model_obs generated values in general
+! Input file can have one qc field, none, or more. note that read_obs_seq_header
+! does NOT return the actual metadata values, which would be helpful in trying
+! to decide if we need to add copies or qcs.
call read_obs_seq_header(obs_sequence_in_name, tnum_copies, tnum_qc, tnum_obs, tmax_num_obs, &
obs_seq_file_id, obs_seq_read_format, pre_I_format, close_the_file = .true.)
-if(tnum_qc == 0) then
- input_qc_index = 1
- DART_qc_index = 2
- ! Need 2 new qc fields, for dummy incoming qc and for the DART qc
- qc_num_inc = 2
- ! original code
- !input_qc_index = 0
- !DART_qc_index = 1
-else if(tnum_qc == 1) then
- input_qc_index = 1
- DART_qc_index = 2
- ! Need 1 new qc field for the DART quality control
- qc_num_inc = 1
+
+! if there are less than 2 incoming qc fields, we will need
+! to make at least 2 (one for the dummy data qc and one for
+! the dart qc).
+if (tnum_qc < 2) then
+ qc_num_inc = 2 - tnum_qc
else
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list