[Dart-dev] [3937] DART/trunk/obs_sequence: Added a copy routine to the obs_sequence_mod that only copies
nancy at ucar.edu
nancy at ucar.edu
Fri Jun 19 12:59:28 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090619/3a30df1a/attachment-0001.html
-------------- next part --------------
Modified: DART/trunk/obs_sequence/obs_sequence_mod.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-06-19 18:53:19 UTC (rev 3936)
+++ DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-06-19 18:59:27 UTC (rev 3937)
@@ -66,7 +66,7 @@
public :: obs_type, init_obs, destroy_obs, get_obs_def, set_obs_def, &
get_obs_values, set_obs_values, replace_obs_values, get_qc, set_qc, &
read_obs, write_obs, replace_qc, interactive_obs, copy_obs, assignment(=), &
- get_obs_key
+ get_obs_key, copy_partial_obs
! Public interfaces for obs covariance modeling
public :: obs_cov_type
@@ -1102,32 +1102,43 @@
type(obs_sequence_type), intent(in) :: seq
character(len = *), intent(in) :: file_name
-integer :: i, file_id
+integer :: i, file_id, rc
integer :: have(max_obs_kinds)
+character(len=11) :: useform
-! Figure out which of the total possible kinds (really types) exist in this
-! sequence, and set the array values to 0 for no, 1 for yes.
-call set_used_kinds(seq, have)
+if(write_binary_obs_sequence) then
+ useform = 'unformatted'
+else
+ useform = 'formatted'
+endif
-! Open the file
+! Open the file. nsc - why is this not using open_file()?
file_id = get_unit()
+write(msg_string, *) 'opening '// trim(useform) // ' file ',trim(file_name)
+call error_handler(E_MSG,'write_obs_seq',msg_string)
+
+open(unit = file_id, file = file_name, form = useform, &
+ action='write', position='rewind', iostat=rc)
+if (rc /= 0) then
+ write(msg_string, *) 'unable to create file '//trim(file_name)
+ call error_handler(E_ERR,'write_obs_seq',msg_string,source,revision,revdate)
+endif
+
+! Write the initial string for help in figuring out binary
if(write_binary_obs_sequence) then
- write(msg_string, *) 'opening unformatted file ',trim(file_name)
- call error_handler(E_MSG,'write_obs_seq',msg_string,source,revision,revdate)
- open(unit = file_id, file = file_name, form = "unformatted")
- ! Write the initial string for help in figuring out binary
write(file_id) 'obs_sequence'
- call write_obs_kind(file_id, 'unformatted', have)
else
- write(msg_string, *) 'opening formatted file ',trim(file_name)
- call error_handler(E_MSG,'write_obs_seq',msg_string,source,revision,revdate)
- open(unit = file_id, file = file_name)
- ! Write the initial string for help in figuring out binary
write(file_id, *) 'obs_sequence'
- call write_obs_kind(file_id, use_list=have)
endif
+! Figure out which of the total possible kinds (really types) exist in this
+! sequence, and set the array values to 0 for no, 1 for yes.
+call set_used_kinds(seq, have)
+
+! Write the TOC, with only the kinds that exist in this seq.
+call write_obs_kind(file_id, useform, have)
+
! First inefficient ugly pass at writing an obs sequence, need to
! update for storage size. CHANGE - use num_obs for the max_num_obs, to
! limit the amount of memory needed when this sequence is read in.
@@ -1142,7 +1153,7 @@
if(write_binary_obs_sequence) then
write(file_id) seq%copy_meta_data(i)
else
- write(file_id, '(a129)') seq%copy_meta_data(i)
+ write(file_id, '(a)') seq%copy_meta_data(i)
endif
end do
@@ -1150,7 +1161,7 @@
if(write_binary_obs_sequence) then
write(file_id) seq%qc_meta_data(i)
else
- write(file_id, '(a129)') seq%qc_meta_data(i)
+ write(file_id, '(a)') seq%qc_meta_data(i)
endif
end do
@@ -1169,7 +1180,7 @@
call close_file(file_id)
write(msg_string, *) 'closed file '//trim(file_name)
-call error_handler(E_MSG,'write_obs_seq',msg_string,source,revision,revdate)
+call error_handler(E_MSG,'write_obs_seq',msg_string)
end subroutine write_obs_seq
@@ -2178,6 +2189,59 @@
end subroutine destroy_obs
+!-----------------------------------------------------
+
+subroutine copy_partial_obs(obs1, obs2, numcopies, copylist, &
+ numqc, qclist)
+
+! Copy from obs2 to obs1, the entire contents of the
+! obs def, but only the copies and qcs listed (in order)
+
+type(obs_type), intent(inout) :: obs1
+type(obs_type), intent(in) :: obs2
+integer, intent(in) :: numcopies, copylist(:), numqc, qclist(:)
+
+integer :: i
+
+! this needs idiotproofing - e.g. bad indices in the copylist
+! or qc list -- without too much expense in time.
+
+obs1%key = obs2%key
+call copy_obs_def(obs1%def, obs2%def)
+
+if (associated(obs1%values)) then
+ if (size(obs1%values) /= numcopies) then
+ deallocate(obs1%values)
+ allocate(obs1%values(numcopies))
+ endif
+else
+ allocate(obs1%values(numcopies))
+endif
+
+if (associated(obs1%qc)) then
+ if (size(obs1%qc) /= numqc) then
+ deallocate(obs1%qc)
+ allocate(obs1%qc(numqc))
+ endif
+else
+ allocate(obs1%qc(numqc))
+endif
+
+do i = 1, numcopies
+ ! error checks here? if copylist(i) > numcopies or < 1, err out
+ obs1%values(i) = obs2%values(copylist(i))
+enddo
+do i = 1, numqc
+ ! and here? if qclist(i) > numqc or < 1, err out
+ obs1%qc(i) = obs2%qc(qclist(i))
+enddo
+
+obs1%prev_time = obs2%prev_time
+obs1%next_time = obs2%next_time
+obs1%cov_group = obs2%cov_group
+
+end subroutine copy_partial_obs
+
!-------------------------------------------------
subroutine get_obs_def(obs, obs_def)
Modified: DART/trunk/obs_sequence/obs_sequence_tool.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-06-19 18:53:19 UTC (rev 3936)
+++ DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-06-19 18:59:27 UTC (rev 3937)
@@ -13,11 +13,11 @@
! $Revision
! $Date
-use types_mod, only : r8, missing_r8
+use types_mod, only : r8, missing_r8, metadatalength, obstypelength
use utilities_mod, only : timestamp, register_module, initialize_utilities, &
find_namelist_in_file, check_namelist_read, &
error_handler, E_ERR, E_MSG, nmlfileunit, &
- do_nml_file, do_nml_term
+ do_nml_file, do_nml_term, get_next_filename
use location_mod, only : location_type, get_location, set_location2, &
LocationName !! , vert_is_height
use obs_def_mod, only : obs_def_type, get_obs_def_time, get_obs_kind, &
@@ -36,7 +36,7 @@
destroy_obs, destroy_obs_sequence, &
delete_seq_head, delete_seq_tail, &
get_num_key_range, delete_obs_by_typelist, &
- get_obs_key, &
+ get_obs_key, copy_partial_obs, &
delete_obs_from_seq, get_next_obs_from_key, &
delete_obs_by_qc, delete_obs_by_copy, &
select_obs_by_location
@@ -49,7 +49,8 @@
revdate = "$Date$"
type(obs_sequence_type) :: seq_in, seq_out
-type(obs_type) :: obs, prev_obs, next_obs, new_obs
+type(obs_type) :: obs_in, next_obs_in
+type(obs_type) :: obs_out, prev_obs_out
logical :: is_there_one, is_this_last
integer :: size_seq_in, num_copies_in, num_qc_in
integer :: size_seq_out, num_copies_out, num_qc_out
@@ -68,19 +69,23 @@
! Namelist input with default values
! max_num_input_files : maximum number of input sequence files to be processed
-integer, parameter :: max_num_input_files = 50
-integer :: num_input_files = 1
+integer, parameter :: max_num_input_files = 500
+integer :: num_input_files = 0
! lazy, pick a big number
integer, parameter :: max_obs_input_types = 500
-character(len = 32) :: obs_types(max_obs_input_types)
+character(len = obstypelength) :: obs_types(max_obs_input_types) = ''
logical :: restrict_by_obs_type
integer :: num_obs_input_types
logical :: restrict_by_location
logical :: restrict_by_qc, restrict_by_copy, restrict_by_height
+logical :: editing_obs, matching_copy_metadata, matching_qc_metadata
+integer :: matching_copy_limit = 0
+integer :: matching_qc_limit = 0
-character(len = 129) :: filename_seq(max_num_input_files) = 'obs_seq.out'
+character(len = 129) :: filename_seq(max_num_input_files) = ''
+character(len = 129) :: filename_seq_list = ''
character(len = 129) :: filename_out = 'obs_seq.processed'
logical :: process_file(max_num_input_files)
@@ -92,7 +97,8 @@
integer :: last_obs_seconds = -1
! must be location independent...
type(location_type) :: min_loc, max_loc
-real(r8) :: min_box(4) = missing_r8, max_box(4) = missing_r8
+real(r8) :: min_box(4) = missing_r8
+real(r8) :: max_box(4) = missing_r8
real(r8) :: min_lat = -90.0_r8
real(r8) :: max_lat = 90.0_r8
real(r8) :: min_lon = 0.0_r8
@@ -102,21 +108,38 @@
real(r8) :: max_qc = missing_r8
real(r8) :: min_copy = missing_r8
real(r8) :: max_copy = missing_r8
-character(len = 32) :: copy_type = ''
-character(len=129) :: qc_metadata = ''
-character(len=129) :: copy_metadata = ''
+character(len = obstypelength) :: copy_type = ''
+character(len=metadatalength) :: copy_metadata = ''
+character(len=metadatalength) :: qc_metadata = ''
+! 256 is an arb max of number of copies for data and qc
+logical :: edit_copy_metadata = .false.
+character(len=metadatalength) :: new_copy_metadata(256) = ''
+logical :: edit_copies = .false.
+integer :: new_copy_index(256) = -1
+integer :: copy_index_len = 0
+logical :: edit_qc_metadata = .false.
+character(len=metadatalength) :: new_qc_metadata(256) = ''
+logical :: edit_qcs = .false.
+integer :: new_qc_index(256) = -1
+integer :: qc_index_len = 0
+character(len=metadatalength) :: synonymous_copy_list(256) = ''
+character(len=metadatalength) :: synonymous_qc_list(256) = ''
logical :: keep_types = .true.
logical :: print_only = .false.
logical :: gregorian_cal = .true.
real(r8) :: min_gps_height = missing_r8
-namelist /obs_sequence_tool_nml/ num_input_files, filename_seq, filename_out, &
+namelist /obs_sequence_tool_nml/ &
+ num_input_files, filename_seq, filename_seq_list, filename_out, &
first_obs_days, first_obs_seconds, last_obs_days, last_obs_seconds, &
- obs_types, keep_types, min_box, max_box, print_only, &
- min_lat, max_lat, min_lon, max_lon, min_qc, max_qc, qc_metadata, &
- min_copy, max_copy, copy_metadata, copy_type, gregorian_cal, &
- min_gps_height
+ obs_types, keep_types, min_box, max_box, print_only, &
+ min_lat, max_lat, min_lon, max_lon, min_qc, max_qc, qc_metadata, &
+ min_copy, max_copy, copy_metadata, copy_type, gregorian_cal, &
+ min_gps_height, edit_copy_metadata, new_copy_index, &
+ edit_qc_metadata, new_qc_index, synonymous_copy_list, &
+ synonymous_qc_list, edit_copies, edit_qcs, new_copy_metadata, &
+ new_qc_metadata
!----------------------------------------------------------------
! Start of the program:
@@ -127,29 +150,30 @@
call obs_seq_modules_used()
-! Initialize input obs_seq filenames and obs_types before reading namelist.
-do i = 1, max_num_input_files
- filename_seq(i) = ""
-enddo
-do i = 1, max_obs_input_types
- obs_types(i) = ""
-enddo
-
! Read the namelist entry
call find_namelist_in_file("input.nml", "obs_sequence_tool_nml", iunit)
read(iunit, nml = obs_sequence_tool_nml, iostat = io)
call check_namelist_read(iunit, io, "obs_sequence_tool_nml")
-if(num_input_files .gt. max_num_input_files) then
- call error_handler(E_ERR,'obs_sequence_tool', &
- 'num_input_files > max_num_input_files. change max_num_input_files in source file', &
- source,revision,revdate)
-endif
-
! Record the namelist values used for the run ...
if (do_nml_file()) write(nmlfileunit, nml=obs_sequence_tool_nml)
if (do_nml_term()) write( * , nml=obs_sequence_tool_nml)
+! ok, here's the new logic:
+! if the user specifies neither filename_seq nor filename_seq_list, we
+! default to trying 'obs_seq.out' and num files = 1.
+! if the user specifies both, it's an error.
+! if the user gives a filelist, we make sure the length is not more
+! than maxfiles and read it into the explicit list and continue.
+! if num_input_files = 0, we count up the non-null items in the list
+! and set num_input_files to that.
+! if the user specifies num_input_files but it doesn't match the list length,
+! we give an error (and maybe suggest 0 if they don't want to keep
+! updating the num_input_files.)
+
+call handle_filenames(filename_seq, filename_seq_list, num_input_files)
+
+
! if you are not using a gregorian cal, set this to false in the namelist.
! if users need it, we could add a calendar type integer to the namelist,
! if users want to specify a particular calendar which is not gregorian.
@@ -300,12 +324,24 @@
restrict_by_copy = .false.
endif
-! Read header information for the sequences to see if we need
-! to accomodate additional copies or qc values from subsequent sequences.
-! Also, calculate how many observations to be added to the first sequence.
-num_copies_out = 0
-num_qc_out = 0
-size_seq_out = 0
+! See if the user is requesting edit of either the number of copies
+! or qc values, or reordering them.
+editing_obs = .false.
+if (edit_copies) then
+ ! count up how many copies and qc there are going out.
+ do i = 1, size(new_copy_index)
+ if (new_copy_index(i) < 0) exit
+ copy_index_len = copy_index_len + 1
+ enddo
+ editing_obs = .true.
+endif
+if (edit_qcs) then
+ do i = 1, size(new_qc_index)
+ if (new_qc_index(i) < 0) exit
+ qc_index_len = qc_index_len + 1
+ enddo
+ editing_obs = .true.
+endif
! check to see if we are going to trim the sequence by time
if(first_obs_seconds >= 0 .or. first_obs_days >= 0) then
@@ -327,6 +363,36 @@
endif
endif
+! check to see if user gave us lists of metadata that should be
+! considered the same across different obs_seq files. count the
+! length of each list. used in the compare_metadata() routine.
+matching_copy_metadata = .false.
+if (synonymous_copy_list(1) /= '') then
+ matching_copy_metadata = .true.
+ do i = 1, size(synonymous_copy_list)
+ if (synonymous_copy_list(i) == '') exit
+ matching_copy_limit = i
+ enddo
+endif
+matching_qc_metadata = .false.
+if (synonymous_qc_list(1) /= '') then
+ matching_qc_metadata = .true.
+ do i = 1, size(synonymous_qc_list)
+ if (synonymous_qc_list(i) == '') exit
+ matching_qc_limit = i
+ enddo
+endif
+
+
+! end of namelist processing and setup
+
+! Read header information for the sequences to see if we need
+! to accomodate additional copies or qc values from subsequent sequences.
+! Also, calculate how many observations to be added to the first sequence.
+num_copies_out = 0
+num_qc_out = 0
+size_seq_out = 0
+
! TWO PASS algorithm; open each file, trim it if requested, and count
! the number of actual observations. then the output file can be
! created with the correct size, and as observations are put into it
@@ -368,10 +434,55 @@
size_seq_in = remaining_obs_count
endif
+ ! First sequence with valid data gets to set up the format of the output.
+ ! This only matters if we keep the 'non-exact match of metadata' in the
+ ! compare_metadata() subroutine below.
if ( first_seq < 0 ) then
first_seq = i
- num_copies_out = num_copies_in
- num_qc_out = num_qc_in
+ ! set up allowing number of copies or qcs to be changed or reordered.
+ if (edit_copies) then
+ ! do a tiny bit of sanity checking here.
+ do j = 1, copy_index_len
+ if ((new_copy_index(j) > num_copies_in) .or. &
+ (new_copy_index(j) < 1)) then
+ write(msgstring,*)'new_copy_index values must be between 1 and ', num_copies_in
+ call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
+ source,revision,revdate)
+
+ endif
+ enddo
+ if (copy_index_len > num_copies_in) then
+ write(msgstring,*)'WARNING: more output copies than input; some are being replicated'
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
+ endif
+ num_copies_out = copy_index_len
+ else
+ num_copies_out = num_copies_in
+ do j = 1, num_copies_in
+ new_copy_index(j) = j
+ enddo
+ endif
+ if (edit_qcs) then
+ do j = 1, qc_index_len
+ if ((new_qc_index(j) > num_qc_in) .or. &
+ (new_qc_index(j) < 1)) then
+ write(msgstring,*)'new_qc_index values must be between 1 and ', num_qc_in
+ call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
+ source,revision,revdate)
+
+ endif
+ enddo
+ if (qc_index_len > num_qc_in) then
+ write(msgstring,*)'WARNING: more output qcs than input; some are being replicated'
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
+ endif
+ num_qc_out = qc_index_len
+ else
+ num_qc_out = num_qc_in
+ do j = 1, num_qc_in
+ new_qc_index(j) = j
+ enddo
+ endif
size_seq_out = size_seq_in
else
size_seq_out = size_seq_out + size_seq_in
@@ -392,10 +503,10 @@
call error_handler(E_MSG,' ',' ')
! Initialize individual observation variables
-call init_obs( obs, num_copies_out, num_qc_out)
-call init_obs( new_obs, num_copies_out, num_qc_out)
-call init_obs(next_obs, num_copies_out, num_qc_out)
-call init_obs(prev_obs, num_copies_out, num_qc_out)
+call init_obs( obs_in, num_copies_in, num_qc_in )
+call init_obs(next_obs_in, num_copies_in, num_qc_in )
+call init_obs( obs_out, num_copies_out, num_qc_out)
+call init_obs(prev_obs_out, num_copies_out, num_qc_out)
total_num_inserted = 0
@@ -426,11 +537,49 @@
if (first_seq < 0) then
call init_obs_sequence(seq_out, num_copies_out, num_qc_out, size_seq_out)
do j=1, num_copies_out
- meta_data = get_copy_meta_data(seq_in, j)
+ if (edit_copy_metadata) then
+ meta_data = new_copy_metadata(j)
+ write(msgstring, *)'replacing original copy meta_data: ' // &
+ trim(get_copy_meta_data(seq_in, new_copy_index(j)))
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ write(msgstring, *)' with: ' // trim(meta_data)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ if (new_copy_index(j) /= j) then
+ write(msgstring, *)'original index was ', new_copy_index(j), ', now ', j
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ endif
+ else
+ meta_data = get_copy_meta_data(seq_in, new_copy_index(j))
+ if (new_copy_index(j) /= j) then
+ write(msgstring, *)'copy with meta_data: ' // trim(meta_data)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ write(msgstring, *)'moved from index ', new_copy_index(j), ' to ', j
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ endif
+ endif
call set_copy_meta_data(seq_out, j, meta_data)
enddo
do j=1, num_qc_out
- meta_data = get_qc_meta_data(seq_in, j)
+ if (edit_qc_metadata) then
+ meta_data = new_qc_metadata(j)
+ write(msgstring, *)'replacing original qc meta_data: ' // &
+ trim(get_qc_meta_data(seq_in, new_qc_index(j)))
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ write(msgstring, *)' with: ' // trim(meta_data)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ if (new_qc_index(j) /= j) then
+ write(msgstring, *)'original index was ', new_qc_index(j), ' now ', j
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ endif
+ else
+ meta_data = get_qc_meta_data(seq_in, new_qc_index(j))
+ if (new_qc_index(j) /= j) then
+ write(msgstring, *)'qc with meta_data: ' // trim(meta_data)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ write(msgstring, *)'moved from index ', new_qc_index(j), ' now ', j
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
+ endif
+ endif
call set_qc_meta_data(seq_out, j, meta_data)
enddo
first_seq = i
@@ -458,20 +607,27 @@
! Must pass a copy of incoming obs to insert_obs_in_seq.
!--------------------------------------------------------------
num_inserted = 0
- is_there_one = get_first_obs(seq_in, obs)
+ is_there_one = get_first_obs(seq_in, obs_in)
if ( is_there_one ) then
- new_obs = obs ! obs records position in seq_out
+ if (editing_obs) then
+ ! copy subset or reordering
+ call copy_partial_obs(obs_out, obs_in, &
+ num_copies_out, new_copy_index, &
+ num_qc_out, new_qc_index)
+ else
+ obs_out = obs_in
+ endif
-!#! call change_variance(new_obs)
+!#! call change_variance(obs_out)
- call insert_obs_in_seq(seq_out, new_obs) ! new_obs linked list info changes
+ call insert_obs_in_seq(seq_out, obs_out) ! new_obs linked list info changes
- prev_obs = new_obs ! records new position in seq_in
+ prev_obs_out = obs_out ! records new position in seq_out
num_inserted = num_inserted + 1
- call get_next_obs(seq_in, obs, next_obs, is_this_last)
+ call get_next_obs(seq_in, obs_in, next_obs_in, is_this_last)
ObsLoop : do while ( .not. is_this_last)
if (print_every > 0) then
@@ -480,22 +636,30 @@
endif
endif
- obs = next_obs ! essentially records position in seq_out
- new_obs = obs ! will be modified w/ position in seq_in
+ obs_in = next_obs_in ! essentially records position in seq_out
-!#! call change_variance(new_obs)
+ if (editing_obs) then
+ ! copy subset or reordering
+ call copy_partial_obs(obs_out, obs_in, &
+ num_copies_out, new_copy_index, &
+ num_qc_out, new_qc_index)
+ else
+ obs_out = obs_in
+ endif
+!#! call change_variance(obs_out)
+
! Since the stride through the observation sequence file is always
! guaranteed to be in temporally-ascending order, we can use the
! 'previous' observation as the starting point to search for the
! correct insertion point. This speeds up the insert code a lot.
- call insert_obs_in_seq(seq_out, new_obs, prev_obs)
+ call insert_obs_in_seq(seq_out, obs_out, prev_obs_out)
- prev_obs = new_obs ! update position in seq_in for next insert
+ prev_obs_out = obs_out ! update position in seq_in for next insert
num_inserted = num_inserted + 1
- call get_next_obs(seq_in, obs, next_obs, is_this_last)
+ call get_next_obs(seq_in, obs_in, next_obs_in, is_this_last)
enddo ObsLoop
@@ -539,10 +703,10 @@
! Time to clean up
call destroy_obs_sequence(seq_out)
-call destroy_obs( obs)
-call destroy_obs( new_obs)
-call destroy_obs(next_obs)
-call destroy_obs(prev_obs)
+call destroy_obs( obs_in )
+call destroy_obs(next_obs_in )
+call destroy_obs( obs_out)
+call destroy_obs(prev_obs_out)
call timestamp(source,revision,revdate,'end')
@@ -571,25 +735,36 @@
!
! The messages might be a bit misleading 'warning', 'warning', 'error' ...
!
+! FIXME: this routine uses several globals now from the namelist that
+! should be arguments to this routine. i'm being lazy (or expedient) here
+! but this should be fixed at some point soon.
+!
+! and now there is an additional restriction -- if editing the copies or
+! qc entries, seq1 has already been edited and only 2 needs the editing
+! applied. before they were completely symmetric.
type(obs_sequence_type), intent(IN) :: seq1, seq2
character(len=*), optional :: fname1, fname2
integer :: num_copies1, num_qc1
integer :: num_copies2, num_qc2
-integer :: num_copies , num_qc, i
+integer :: num_copies , num_qc, i, j
+logical :: have_match1, have_match2
character(len=129) :: str1, str2
character(len=255) :: msgstring1, msgstring2
+num_copies1 = get_num_copies(seq1)
num_qc1 = get_num_qc( seq1)
-num_copies1 = get_num_copies(seq1)
+num_copies2 = get_num_copies(seq2)
num_qc2 = get_num_qc( seq2)
-num_copies2 = get_num_copies(seq2)
-num_copies = num_copies1
-num_qc = num_qc1
+if (edit_copies) num_copies2 = copy_index_len
+if (edit_qcs ) num_qc2 = qc_index_len
+num_copies = num_copies2
+num_qc = num_qc2
+
! get this ready in case we have to use it
if (present(fname1) .and. present(fname2)) then
write(msgstring1,*)'Sequence files ', trim(fname1), ' and ', trim(fname2), &
@@ -614,36 +789,157 @@
call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
endif
-MetaDataLoop : do i=1, num_copies
- str1 = trim(adjustl(get_copy_meta_data(seq1,i)))
- str2 = trim(adjustl(get_copy_meta_data(seq2,i)))
+! watch the code flow in this loop and the one below it.
+! the smoothest code order is to determine what the strings are,
+! and then try different things to match them. if a match is found,
+! cycle. if you get to the bottom of the loop, there is no match
+! and a single set of (fatal) error messages is called.
+CopyMetaData : do i=1, num_copies
+ str1 = get_copy_meta_data(seq1,i)
+ if (edit_copy_metadata) then
+ str2 = new_copy_metadata(i)
+ else
+ str2 = get_copy_meta_data(seq2, new_copy_index(i))
+ endif
+
+ ! easy case - they match. cycle to next copy.
if( str1 == str2 ) then
- write(msgstring2,*)'metadata ',trim(adjustl(str1)), ' in both.'
+ write(msgstring2,*)'metadata ',trim(str1), ' in both.'
call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- else
- write(msgstring2,*)'metadata value mismatch. seq1: ', trim(adjustl(str1))
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- write(msgstring2,*)'metadata value mismatch. seq2: ', trim(adjustl(str2))
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+ cycle CopyMetaData
endif
-enddo MetaDataLoop
+ ! see if user provided a list of metadata strings that are
+ ! the same values and can be considered a match.
+ if (matching_copy_metadata) then
+ have_match1 = .false.
+ do j=1, matching_copy_limit
+ if (trim(str1) == trim(synonymous_copy_list(j))) then
+ have_match1 = .true.
+ exit
+ endif
+ enddo
+
+ have_match2 = .false.
+ do j=1, matching_copy_limit
+ if (trim(str2) == trim(synonymous_copy_list(j))) then
+ have_match2 = .true.
+ exit
+ endif
+ enddo
+
+ ! if both are true, you found both strings in the list and
+ ! it is ok to proceed.
+ if (have_match1 .and. have_match2) then
+ write(msgstring2,*)'different copy metadata strings ok because both on synonymous list'
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'one is: ', trim(str1)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'one is: ', trim(str2)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ cycle CopyMetaData
+ endif
+
+ ! if no match, fall out of the if.
+ endif
+
+ !! FIXME: this could be dangerous - it allows any metadata string with
+ !! the substring 'observation' to match any other. if there are multiple
+ !! strings with 'observation', it will allow them to match. for now it
+ !! allows 'NCEP BUFR observations' to match 'observation', for example,
+ !! but it's dangerous.
+ !if ((index(str1, 'observation') > 0) .and. &
+ ! (index(str2, 'observation') > 0)) then
+ ! write(msgstring2,*)'observation metadata in both ',trim(str1), ' and ', trim(str2)
+ ! call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! write(msgstring2,*)'ALLOWING NON-EXACT MATCH'
+ ! call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! cycle CopyMetaData
+ !! end FIXME
+
+ ! if you get here, the metadata is not the same and the user has not
+ ! given us strings that are ok to match. fail.
+ write(msgstring2,*)'metadata value mismatch. seq1: ', trim(str1)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'metadata value mismatch. seq2: ', trim(str2)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+
+enddo CopyMetaData
+
QCMetaData : do i=1, num_qc
- str1 = trim(adjustl(get_qc_meta_data(seq1,i)))
- str2 = trim(adjustl(get_qc_meta_data(seq2,i)))
+ str1 = get_qc_meta_data(seq1,i)
+ if (edit_qc_metadata) then
+ str2 = new_qc_metadata(i)
+ else
+ str2 = get_qc_meta_data(seq2, new_qc_index(i))
+ endif
+
+
+ ! easy case - they match. cycle to next copy.
if( str1 == str2 ) then
- write(msgstring2,*)'qc metadata ', trim(adjustl(str1)), ' in both.'
+ write(msgstring2,*)'metadata ',trim(str1), ' in both.'
call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- else
- write(msgstring2,*)'qc metadata value mismatch. seq1: ', trim(adjustl(str1))
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- write(msgstring2,*)'qc metadata value mismatch. seq2: ', trim(adjustl(str2))
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
- call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+ cycle QCMetaData
endif
+
+ ! see if user provided a list of metadata strings that are
+ ! the same values and can be considered a match.
+ if (matching_copy_metadata) then
+ have_match1 = .false.
+ do j=1, matching_qc_limit
+ if (trim(str1) == trim(synonymous_qc_list(j))) then
+ have_match1 = .true.
+ exit
+ endif
+ enddo
+
+ have_match2 = .false.
+ do j=1, matching_qc_limit
+ if (trim(str2) == trim(synonymous_qc_list(j))) then
+ have_match2 = .true.
+ exit
+ endif
+ enddo
+
+ ! if both are true, you found both strings in the list and
+ ! it is ok to proceed.
+ if (have_match1 .and. have_match2) then
+ write(msgstring2,*)'different qc metadata strings ok because both on synonymous list'
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'one is: ', trim(str1)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'one is: ', trim(str2)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ cycle QCMetaData
+ endif
+
+ ! if no match, fall out of the if.
+ endif
+
+ !! FIXME: this is even more dangerous than the obs - better to make the
+ !! user give an explicit list of strings that are ok to match. but here
+ !! is the code if you wanted to make it more mindless.
+ !if (((index(str1, 'QC') > 0).or.(index(str1, 'quality control') > 0)).and. &
+ ! ((index(str2, 'QC') > 0).or.(index(str2, 'quality control') > 0))) then
+ ! write(msgstring2,*)'QC metadata in both ',trim(str1), ' and ', trim(str2)
+ ! call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! write(msgstring2,*)'ALLOWING NON-EXACT MATCH'
+ ! call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ ! cycle QCMetaData
+ !endif
+ !! end FIXME
+
+ ! if you get here, the metadata is not the same and the user has not
+ ! given us strings that are ok to match. fail.
+ write(msgstring2,*)'qc metadata value mismatch. seq1: ', trim(str1)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ write(msgstring2,*)'qc metadata value mismatch. seq2: ', trim(str2)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
+ call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
+
enddo QCMetaData
end subroutine compare_metadata
@@ -670,7 +966,7 @@
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
- ' are before first_obs_days:first_obs_seconds'
+ ' are before first_obs time'
call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
@@ -684,7 +980,7 @@
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
- ' are after last_obs_days:last_obs_seconds'
+ ' are after last_obs time'
call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
@@ -1026,21 +1322,22 @@
num_qc = get_num_qc( seq1)
if ( num_copies < 0 .or. num_qc < 0 ) then
+ write(msgstring1,*)' illegal copy or obs count in file '//trim(fname1)
call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
endif
MetaDataLoop : do i=1, num_copies
- str1 = trim(adjustl(get_copy_meta_data(seq1,i)))
+ str1 = get_copy_meta_data(seq1,i)
- write(msgstring1,*)'Data Metadata: ',trim(adjustl(str1))
+ write(msgstring1,*)'Data Metadata: ',trim(str1)
call error_handler(E_MSG, '', msgstring1)
enddo MetaDataLoop
QCMetaData : do i=1, num_qc
- str1 = trim(adjustl(get_qc_meta_data(seq1,i)))
+ str1 = get_qc_meta_data(seq1,i)
- write(msgstring1,*)' QC Metadata: ', trim(adjustl(str1))
+ write(msgstring1,*)' QC Metadata: ', trim(str1)
call error_handler(E_MSG, '', msgstring1)
enddo QCMetaData
@@ -1170,7 +1467,99 @@
end subroutine select_gps_by_height
+!---------------------------------------------------------------------
+subroutine handle_filenames(filename_seq, filename_seq_list, num_input_files)
+! sort out the input lists, set the length if not given by user,
+! make sure what's specified is consistent.
+character(len=*), intent(inout) :: filename_seq(:)
+character(len=*), intent(in) :: filename_seq_list
+integer, intent(inout) :: num_input_files
+integer :: index
+logical :: from_file
+character(len=32) :: source
+
+! ok, here's the new logic:
+! if the user specifies neither filename_seq nor filename_seq_list, we
+! default to trying 'obs_seq.out' and num files = 1.
+! if the user specifies both, it's an error.
+! if the user gives a filelist, we make sure the length is not more
+! than maxfiles and read it into the explicit list and continue.
+! if num_input_files = 0, we count up the non-null items in the list
+! and set num_input_files to that.
+! if the user specifies num_input_files but it doesn't match the list length,
+! we give an error (and maybe suggest 0 if they don't want to keep
+! updating the num_input_files.)
+
+! default case - input file is 'obs_seq.out' and count is 1.
+if (filename_seq(1) == '' .and. filename_seq_list == '') then
+
+ if (num_input_files /= 0 .and. num_input_files /= 1) then
+ call error_handler(E_ERR,'obs_sequence_tool', &
+ 'if no filenames specified, num_input_files must be 0 or 1', &
+ source,revision,revdate)
+ endif
+
+ num_input_files = 1
+ filename_seq(1) = 'obs_seq.out'
+ return
+endif
+
+! make sure the namelist specifies one or the other but not both
+if (filename_seq(1) /= '' .and. filename_seq_list /= '') then
+ call error_handler(E_ERR,'obs_sequence_tool', &
+ 'cannot specify both filename_seq and filename_seq_list', &
+ source,revision,revdate)
+endif
+
+! if they have specified a file which contains a list, read it into
+! the filename_seq array and set the count.
+if (filename_seq_list /= '') then
+ source = 'filename_seq_list'
+ from_file = .true.
+else
+ source = 'filename_seq'
+ from_file = .false.
+endif
+
+do index = 1, max_num_input_files
+ if (from_file) &
+ filename_seq(index) = get_next_filename(filename_seq_list, index)
+
+ if (filename_seq(index) == '') then
+ if (index == 1) then
+ call error_handler(E_ERR,'obs_sequence_tool', &
+ trim(source)//' contains no filenames', &
+ source,revision,revdate)
+ endif
+ ! leaving num_input_files unspecified (or set to 0) means use
+ ! whatever number of files is in the list.
+ if (num_input_files == 0) then
+ num_input_files = index - 1
+ return
+ else
+ ! if they do give a count, make it match.
+ if (num_input_files == (index - 1)) return
+
+ write(msgstring, *) 'if num_input_files is 0, the number of files will be automatically computed'
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
+ write(msgstring, *) 'if num_input_files is not 0, it must match the number of filenames specified'
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
+ write(msgstring, *) 'num_input_files is ', num_input_files, &
+ ' but '//trim(source)//' has filecount ', index - 1
+ call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
+ source,revision,revdate)
+
+ endif
+ endif
+enddo
+
+write(msgstring, *) 'cannot specify more than ',max_num_input_files,' files'
+call error_handler(E_ERR,'obs_sequence_tool', msgstring, &
+ source,revision,revdate)
+
+end subroutine handle_filenames
+
!---------------------------------------------------------------------
!#! subroutine change_variance(this_obs)
!#!
Modified: DART/trunk/obs_sequence/obs_sequence_tool.html
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_tool.html 2009-06-19 18:53:19 UTC (rev 3936)
+++ DART/trunk/obs_sequence/obs_sequence_tool.html 2009-06-19 18:59:27 UTC (rev 3937)
@@ -32,6 +32,11 @@
-->
<DIV ALIGN=CENTER>
+<A HREF="#OVERVIEW">OVERVIEW</A> /
+<A HREF="#Examples">EXAMPLES</A> /
+<A HREF="#Discussion">DISCUSSION</A> /
+<A HREF="#FAQ">FAQ</A> /
+<A HREF="#Building">BUILDING</A> /
<A HREF="#Namelist">NAMELIST</A> /
<A HREF="#FilesUsed">FILES</A> /
<A HREF="#References">REFERENCES</A> /
@@ -52,65 +57,500 @@
<H2>OVERVIEW</H2>
<P>
-This program processes one or more compatible observation sequences
-into a single observation sequence file, optionally selecting ranges
-by time, observation type, or location box.
-The number of input files, their filenames, the output filename, and
-all the optional selection settings are read from the
-<em class=code>&obs_sequence_tool_nml</em> namelist
-in file <em class=file>input.nml</em>.
-<BR><BR>
+DART observation sequence files are stored in a proprietary format.
+This tool makes it easier to manipulate these files, allowing the user
+to subset or combine one or more files into a single output file.
+</P>
+
+<P>
+The tool has many options
+to select subsets of observations by time, type, data value, and location.
+The tool also allows the contents of observations to be
+changed by subsetting and/or reordering the copies and qc entries.
+Files with equivalent data but with different metadata labels
+(e.g. 'NCEP QC' vs. 'QC') can now be merged as well.
+The tool can be run without creating an output file, only printing
+a summary of the counts of each observation type in the input files,
+and it can be used to convert from binary to ASCII and back.
+</P>
+
+<P>
+The actions of the <em class=code>obs_sequence_tool</em> program
+are controlled by a Fortran namelist, read from a file named
+<em class=file>input.nml</em> in the current directory. A detailed
+description of each namelist item is described in the
+<a href="#Namelist">namelist section</a> of this document.
+The names used in this discussion refer to these namelist items.
+</P>
+
+<P>
+The following section contains examples of common usages for this tool.
+Below that are more details about DART observation sequence files, the
+structure of individual observations, and general background information.
+</P>
+
+
@@ Diff output truncated at 40000 characters. @@
More information about the Dart-dev
mailing list