[Dart-dev] [3220] DART/branches/nancy_work/obs_sequence: My test
branch: further updates, not ready for trunk yet.
nancy at subversion.ucar.edu
nancy at subversion.ucar.edu
Fri Feb 8 16:28:35 MST 2008
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20080208/fbbc9b9b/attachment.html
-------------- next part --------------
Modified: DART/branches/nancy_work/obs_sequence/merge_obs_seq.f90
===================================================================
--- DART/branches/nancy_work/obs_sequence/merge_obs_seq.f90 2008-02-08 00:04:53 UTC (rev 3219)
+++ DART/branches/nancy_work/obs_sequence/merge_obs_seq.f90 2008-02-08 23:28:35 UTC (rev 3220)
@@ -39,9 +39,10 @@
integer :: size_seq_in, num_copies_in, num_qc_in
integer :: size_seq_out, num_copies_out, num_qc_out
integer :: num_inserted, iunit, io, i, j, total_num_inserted
-integer :: max_num_obs, file_id, first_seq_index, remaining_obs_count
+integer :: max_num_obs, file_id, remaining_obs_count
+integer :: first_seq
character(len = 129) :: read_format, meta_data
-logical :: pre_I_format, all_gone, first_seq
+logical :: pre_I_format, all_gone
logical :: trim_first, trim_last
character(len = 129) :: msgstring
@@ -122,7 +123,7 @@
endif
if (trim_first .and. trim_last) then
if (first_obs_time > last_obs_time) then
- call error_handler(E_ERR,'merge_obs_seq', 'first time later than last time', &
+ call error_handler(E_ERR,'merge_obs_seq', 'first time cannot be later than last time', &
source,revision,revdate)
endif
endif
@@ -134,8 +135,7 @@
! pass 1:
-first_seq_index = -1
-first_seq = .true.
+first_seq = -1
do i = 1, num_input_files
if ( len(filename_seq(i)) .eq. 0 .or. filename_seq(i) .eq. "" ) then
@@ -143,54 +143,45 @@
'num_input_files and filename_seq mismatch',source,revision,revdate)
endif
- ! if some of the files are completely removed, you can ignore them.
+ ! count up the number of observations we are going to eventually have.
+ ! if all the observations in a file are not part of the linked list, the
+ ! output number of observations might be much smaller than the total size in
+ ! the header. it is slower, but go ahead and read in the entire sequence
+ ! and count up the real number of obs - trim_seq will do the count even if
+ ! it is not trimming in time. this allows us to create an empty obs_seq
+ ! output file of exactly the right size.
+
call read_obs_seq_header(filename_seq(i), num_copies_in, num_qc_in, size_seq_in, &
max_num_obs, file_id, read_format, pre_I_format, close_the_file = .true.)
- ! if not trimming observations, just read the headers because it is faster.
- if (trim_first .or. trim_last) then
- call read_obs_seq(filename_seq(i), 0, 0, 0, seq_in)
- call trim_seq(seq_in, trim_first, first_obs_time, trim_last, last_obs_time, &
- filename_seq(i), .true., remaining_obs_count)
- call destroy_obs_sequence(seq_in)
- if (remaining_obs_count == 0) then
- process_file(i) = .false.
- cycle
- else
- process_file(i) = .true.
- endif
+
+ call read_obs_seq(filename_seq(i), 0, 0, 0, seq_in)
+ call trim_seq(seq_in, trim_first, first_obs_time, trim_last, last_obs_time, &
+ filename_seq(i), .true., remaining_obs_count)
+ call destroy_obs_sequence(seq_in)
+ if (remaining_obs_count == 0) then
+ process_file(i) = .false.
+ cycle
+ else
+ process_file(i) = .true.
size_seq_in = remaining_obs_count
endif
-
- if ( first_seq ) then
- first_seq_index = i
- first_seq = .false.
+ if ( first_seq < 0 ) then
+ first_seq = i
num_copies_out = num_copies_in
num_qc_out = num_qc_in
size_seq_out = size_seq_in
else
- ! some simple error checking here; first, counts must match.
- ! later we will verify that the metadata contents themselves match.
- if (num_copies_out /= num_copies_in) then
- write(msgstring,51) 'Number of data copies inconsistent, ', &
- num_copies_in, ' does not equal ', num_copies_out
- call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
- endif
- if (num_qc_out /= num_qc_in) then
- write(msgstring,51) 'Number of QC copies inconsistent, ', &
- num_qc_in, ' does not equal ', num_qc_out
- call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
- endif
size_seq_out = size_seq_out + size_seq_in
endif
-51 format (A,I4,A,I4)
enddo
-! no valid obs found?
-if (first_seq_index == -1) then
- msgstring = 'All obs in all files are empty or were deleted by the first/last times'
- call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
+! no valid obs found? if the index value is still negative, we are
+! still waiting to process the first one and never found one.
+if (first_seq < 0 .or. size_seq_out == 0) then
+ msgstring = 'All input files are empty or all obs outside the first/last times'
+ call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
endif
! pass 2:
@@ -204,53 +195,62 @@
total_num_inserted = 0
! Read obs seq to be added, and insert obs from it to the output seq
+first_seq = -1
do i = 1, num_input_files
if (.not. process_file(i)) cycle
+ write(msgstring,*) 'Starting to process input sequence file ', trim(filename_seq(i))
+ call error_handler(E_MSG,'merge_obs_seq',msgstring,source,revision,revdate)
+
call read_obs_seq(filename_seq(i), 0, 0, 0, seq_in)
- ! Need to find first obs with appropriate time, delete all earlier ones
+ ! If you get here, there better be observations in this file which
+ ! are going to be used (the process_file flag wouldn't be set otherwise.)
call trim_seq(seq_in, trim_first, first_obs_time, trim_last, last_obs_time, &
filename_seq(i), .false., remaining_obs_count)
+
+ ! This would be an error at this point.
if(remaining_obs_count == 0) then
- call destroy_obs_sequence(seq_in)
- cycle ! skip this input file and do the next (if there are any)
+ call destroy_obs_sequence(seq_in)
+ write(msgstring, *) 'Internal error trying to process file ', trim(filename_seq(i))
+ call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
endif
-
- ! create the output sequence here
- if (i == first_seq_index) then
- call init_obs_sequence(seq_out, num_copies_out, num_qc_out, size_seq_out)
+
+ ! create the output sequence here based on the first input file
+ 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)
+ meta_data = get_copy_meta_data(seq_in, j)
call set_copy_meta_data(seq_out, j, meta_data)
- enddo
+ enddo
do j=1, num_qc_out
- meta_data = get_qc_meta_data(seq_in, j)
+ meta_data = get_qc_meta_data(seq_in, j)
call set_qc_meta_data(seq_out, j, meta_data)
- enddo
+ enddo
+ first_seq = i
+ else
+ ! we have an existing output sequence already. make sure the next one
+ ! is completely compatible.
+
+ ! Compare metadata between the observation sequences.
+ ! This routine exits if they do not match.
+ call compare_metadata(seq_out, seq_in, filename_seq(first_seq), filename_seq(i))
endif
size_seq_out = get_num_key_range(seq_out) !current size of seq_out
size_seq_in = get_num_key_range(seq_in) !current size of seq_in
- ! Compare metadata between the observation sequences.
- ! If it is compatible, keep going, if not ... it will terminate here.
- call compare_metadata(seq_out, seq_in)
-
- call error_handler(E_MSG,'merge_obs_seq', &
- 'Good news - observation sequence files compatible ...', source,revision,revdate)
-
! Getting the time of the last observation in the first sequence to see
! if we can append instead of insert. Appending is MUCH faster.
is_there_one = get_last_obs(seq_out, obs)
if (.not. is_there_one .and. size_seq_out >= 1) then
- call error_handler(E_ERR,'merge_obs_seq', &
- 'BAD news - first obs_seq is neverending ...', source,revision,revdate)
+ write(msgstring, *) 'Internal error trying to get last obs in ', trim(filename_seq(i))
+ call error_handler(E_ERR,'merge_obs_seq',msgstring,source,revision,revdate)
endif
!-------------------------------------------------------------
- ! Start to insert obs from sequence 2 into sequence 1
+ ! Start to insert obs from sequence_in into sequence_out
!
! NOTE: insert_obs_in_seq CHANGES the obs passed in.
! Must pass a copy of incoming obs to insert_obs_in_seq.
@@ -339,8 +339,7 @@
end subroutine merge_obs_seq_modules_used
- subroutine compare_metadata(seq1, seq2)
-! subroutine compare_metadata(seq1, seq2)
+ subroutine compare_metadata(seq1, seq2, fname1, fname2)
!
! This subroutine compares the metadata for two different observation
! sequences and terminates the program if they are not conformable.
@@ -351,12 +350,13 @@
!
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
character(len=129) :: str1, str2
-character(len=255) :: msgstring
+character(len=255) :: msgstring1, msgstring2
num_qc1 = get_num_qc( seq1)
num_copies1 = get_num_copies(seq1)
@@ -367,21 +367,28 @@
num_copies = num_copies1
num_qc = num_qc1
+! 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), &
+ ' are not compatible'
+else
+ msgstring1 = 'Sequence files cannot be merged because they are not compatible'
+endif
+
if ( num_copies1 /= num_copies2 ) then
- write(msgstring,*)'The obs_sequences have incompatible numbers of copies', &
- num_copies1, num_copies2
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
+ write(msgstring2,*)'Different numbers of data copies found: ', &
+ num_copies1, ' vs ', num_copies2
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
num_copies = -1
endif
if ( num_qc1 /= num_qc2 ) then
- write(msgstring,*)'The obs_sequences have incompatible numbers of QC metadata', &
- num_qc1, num_qc2
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
+ write(msgstring2,*)'Different different numbers of QCs found: ', &
+ num_qc1, ' vs ', num_qc2
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
num_qc = -1
endif
if ( num_copies < 0 .or. num_qc < 0 ) then
- call error_handler(E_ERR, 'merge_obs_seq', &
- 'obs_sequence files incompatible ... stopping.', source, revision, revdate)
+ call error_handler(E_ERR, 'merge_obs_seq', msgstring1, source, revision, revdate)
endif
MetaDataLoop : do i=1, num_copies
@@ -389,15 +396,14 @@
str2 = trim(adjustl(get_copy_meta_data(seq2,i)))
if( str1 == str2 ) then
- write(msgstring,*)'metadata ',trim(adjustl(str1)), ' in both.'
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
+ write(msgstring2,*)'metadata ',trim(adjustl(str1)), ' in both.'
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
else
- write(msgstring,*)'metadata seq1 ', trim(adjustl(str1))
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
- write(msgstring,*)'metadata seq2 ', trim(adjustl(str2))
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
- call error_handler(E_ERR, 'merge_obs_seq', &
- 'obs_sequence files incompatible ... stopping.', source, revision, revdate)
+ write(msgstring2,*)'metadata value mismatch. seq1: ', trim(adjustl(str1))
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
+ write(msgstring2,*)'metadata value mismatch. seq2: ', trim(adjustl(str2))
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
+ call error_handler(E_ERR, 'merge_obs_seq', msgstring1, source, revision, revdate)
endif
enddo MetaDataLoop
@@ -406,15 +412,14 @@
str2 = trim(adjustl(get_qc_meta_data(seq2,i)))
if( str1 == str2 ) then
- write(msgstring,*)'qc metadata ', trim(adjustl(str1)), ' in both.'
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
+ write(msgstring2,*)'qc metadata ', trim(adjustl(str1)), ' in both.'
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
else
- write(msgstring,*)'qc metadata seq1 ', trim(adjustl(str1))
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
- write(msgstring,*)'qc metadata seq2 ', trim(adjustl(str2))
- call error_handler(E_MSG, 'merge_obs_seq', msgstring, source, revision, revdate)
- call error_handler(E_ERR, 'merge_obs_seq', &
- 'obs_sequence files incompatible ... stopping.', source, revision, revdate)
+ write(msgstring2,*)'qc metadata value mismatch. seq1: ', trim(adjustl(str1))
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
+ write(msgstring2,*)'qc metadata value mismatch. seq2: ', trim(adjustl(str2))
+ call error_handler(E_MSG, 'merge_obs_seq', msgstring2, source, revision, revdate)
+ call error_handler(E_ERR, 'merge_obs_seq', msgstring1, source, revision, revdate)
endif
enddo QCMetaData
Modified: DART/branches/nancy_work/obs_sequence/obs_sequence_mod.f90
===================================================================
--- DART/branches/nancy_work/obs_sequence/obs_sequence_mod.f90 2008-02-08 00:04:53 UTC (rev 3219)
+++ DART/branches/nancy_work/obs_sequence/obs_sequence_mod.f90 2008-02-08 23:28:35 UTC (rev 3220)
@@ -1,5 +1,5 @@
! Data Assimilation Research Testbed -- DART
-! Copyright 2004-2007, Data Assimilation Research Section
+! Copyright 2004-2008, Data Assimilation Research Section
! University Corporation for Atmospheric Research
! Licensed under the GPL -- www.gpl.org/licenses/gpl.html
@@ -32,9 +32,9 @@
use time_manager_mod, only : time_type, operator(>), operator(<), operator(>=), &
operator(/=), set_time, operator(-), operator(+), &
operator(==)
-use utilities_mod, only : get_unit, close_file, find_namelist_in_file, check_namelist_read, &
- register_module, error_handler, E_ERR, E_WARN, E_MSG, logfileunit, &
- do_output
+use utilities_mod, only : get_unit, close_file, register_module, error_handler, &
+ find_namelist_in_file, check_namelist_read, &
+ E_ERR, E_WARN, E_MSG, logfileunit, do_output
use obs_kind_mod, only : write_obs_kind, read_obs_kind
@@ -50,11 +50,12 @@
get_num_copies, get_num_qc, get_num_obs, get_max_num_obs, get_copy_meta_data, &
get_qc_meta_data, get_next_obs, get_prev_obs, insert_obs_in_seq, &
delete_obs_from_seq, set_copy_meta_data, set_qc_meta_data, get_first_obs, &
- get_last_obs, add_copies, add_qc, write_obs_seq, read_obs_seq, &
- append_obs_to_seq, get_obs_from_key, get_obs_time_range, set_obs, get_time_range_keys, &
+ get_last_obs, add_copies, add_qc, write_obs_seq, read_obs_seq, set_obs, &
+ append_obs_to_seq, get_obs_from_key, get_obs_time_range, get_time_range_keys, &
get_num_times, get_num_key_range, &
static_init_obs_sequence, destroy_obs_sequence, read_obs_seq_header, &
- get_expected_obs, delete_seq_head, delete_seq_tail, get_next_obs_from_key, get_prev_obs_from_key
+ get_expected_obs, delete_seq_head, delete_seq_tail, &
+ get_next_obs_from_key, get_prev_obs_from_key, copy_obs_seq
! Public interfaces for obs
public :: obs_type, init_obs, destroy_obs, get_obs_def, set_obs_def, &
@@ -76,12 +77,12 @@
integer :: num_qc
integer :: num_obs
integer :: max_num_obs
- character(len = 129), pointer :: copy_meta_data(:)
- character(len = 129), pointer :: qc_meta_data(:)
+ character(len = 129), pointer :: copy_meta_data(:) => NULL()
+ character(len = 129), pointer :: qc_meta_data(:) => NULL()
integer :: first_time
integer :: last_time
! integer :: first_avail_time, last_avail_time
- type(obs_type), pointer :: obs(:)
+ type(obs_type), pointer :: obs(:) => NULL()
! What to do about groups
end type obs_sequence_type
@@ -91,8 +92,8 @@
! Do I want to enforce the identity of the particular obs_sequence?
integer :: key
type(obs_def_type) :: def
- real(r8), pointer :: values(:)
- real(r8), pointer :: qc(:)
+ real(r8), pointer :: values(:) => NULL()
+ real(r8), pointer :: qc(:) => NULL()
! Put sort indices directly into the data structure
integer :: prev_time, next_time
integer :: cov_group
@@ -147,7 +148,7 @@
!--------------------------------------------------------------
!WHAT ABOUT PASS THROUGHS TO THE OBS_DEF???
-! WhAT ABOUT copy_obs_sequence similar to read.
+! WHAT ABOUT copy_obs_sequence similar to read.
!-------------------------------------------------
subroutine init_obs_sequence(seq, num_copies, num_qc, &
expected_max_num_obs)
@@ -206,16 +207,28 @@
if ( seq%max_num_obs > 0 ) then
- if (associated(seq%copy_meta_data)) deallocate(seq%copy_meta_data)
- if (associated(seq%qc_meta_data)) deallocate(seq%qc_meta_data)
+ if (associated(seq%copy_meta_data)) then
+ deallocate(seq%copy_meta_data)
+ nullify(seq%copy_meta_data)
+ endif
+ if (associated(seq%qc_meta_data)) then
+ deallocate(seq%qc_meta_data)
+ nullify(seq%qc_meta_data)
+ endif
do i = 1, seq%max_num_obs
+ ! seq%obs is a derived type, not a pointer.
! if (associated(seq%obs(i))) call destroy_obs( seq%obs(i) )
- call destroy_obs( seq%obs(i) )
+ call destroy_obs( seq%obs(i) )
end do
! Also free up the obs storage in the sequence
- if(associated(seq%obs)) deallocate(seq%obs)
+ if(associated(seq%obs)) then
+ deallocate(seq%obs)
+ nullify(seq%obs)
+ else
+ print *, 'destroy_obs_sequence called but seq%obs not associated'
+ endif
seq%first_time = -1
seq%last_time = -1
@@ -313,6 +326,7 @@
call init_obs(obs, 0, 0)
do i = 1, num_obs
+!write(*,*) 'get_expected_obs: i, num_obs = ', i, num_obs
call get_obs_from_key(seq, keys(i), obs)
call get_obs_def(obs, obs_def)
location = get_obs_def_location(obs_def)
@@ -1481,7 +1495,76 @@
end subroutine delete_seq_tail
+!------------------------------------------------------------------
+! Follow the linked list entries to copy only the linked observations
+! from one sequence to the other. TODO: need a way to get the count
+! of the linked observations only.
+
+subroutine copy_obs_seq(oldseq, newseq, time1, time2)
+
+type(obs_sequence_type), intent(in) :: oldseq
+type(obs_sequence_type), intent(out) :: newseq
+type(time_type), optional, intent(in) :: time1, time2
+
+integer :: i, num_copies, num_qc, max_num_obs
+integer :: num_obs, num_real_obs, num_keys, key_bounds(2)
+integer, pointer :: keylist(:)
+type(obs_type) :: obs
+type(time_type) :: first_time, last_time
+logical :: out_of_range
+
+! Get existing header info
+num_copies = get_num_copies(oldseq)
+num_qc = get_num_qc(oldseq)
+num_obs = get_num_obs(oldseq)
+max_num_obs = get_max_num_obs(oldseq)
+
+
+! Really count how many obs are in the linked list, with
+! optional time starts and ends.
+if (present(time1)) then
+ first_time = time1
+else
+ call get_obs_from_key(oldseq, oldseq%first_time, obs)
+ first_time = get_obs_def_time(obs%def)
+endif
+if (present(time2)) then
+ last_time = time2
+else
+ call get_obs_from_key(oldseq, oldseq%last_time, obs)
+ first_time = get_obs_def_time(obs%def)
+endif
+
+call get_obs_time_range(oldseq, first_time, last_time, &
+ key_bounds, num_keys, out_of_range)
+if (out_of_range) then
+ write(msg_string, *) 'All keys out of range'
+ call error_handler(E_ERR, 'copy_obs_seq', msg_string, &
+ source, revision, revdate)
+endif
+
+
+call init_obs_sequence(newseq, num_copies, num_qc, num_keys)
+
+allocate(keylist(num_keys))
+call get_time_range_keys(oldseq, key_bounds, num_keys, keylist)
+
+call init_obs(obs, num_copies, num_qc)
+
+do i=1, num_keys
+ call get_obs_from_key(oldseq, keylist(i), obs)
+ call set_obs(newseq, obs, i)
+enddo
+
+
+! Release the temp storage
+deallocate(keylist)
+call destroy_obs(obs)
+
+end subroutine copy_obs_seq
+
+
!=================================================
! Functions for the obs_type
@@ -1509,19 +1592,21 @@
! To be overloaded with =
-type(obs_type), intent(out) :: obs1
+!type(obs_type), intent(out) :: obs1
+type(obs_type), intent(inout) :: obs1
type(obs_type), intent(in) :: obs2
obs1%key = obs2%key
call copy_obs_def(obs1%def, obs2%def)
-!write(*, *) 'in copy obs'
-!write(*, *) 'size of obs1, obs2 ', size(obs1%values), size(obs2%values)
if (.not.associated(obs1%values) .or. .not.associated(obs1%qc) .or. &
size(obs1%values) /= size(obs2%values) .or. size(obs1%qc) /= size(obs2%qc)) then
if (associated(obs1%values)) deallocate(obs1%values)
if (associated(obs1%qc)) deallocate(obs1%qc)
- !write(*, *) 'allocating in copy_obs'
+ !write(*, *) 'copy_obs: assoc obs1 values, qc = ', associated(obs1%values), associated(obs1%qc)
+ !write(*, *) 'copy_obs: size of obs1, obs2 values = ', size(obs1%values), size(obs2%values)
+ !write(*, *) 'copy_obs: size of obs1, obs2 qc = ', size(obs1%qc), size(obs2%qc)
+ !write(*, *) 'copy_obs: allocating total bytes = ', size(obs2%values) + size(obs2%qc)
allocate(obs1%values(size(obs2%values)), obs1%qc(size(obs2%qc)))
endif
obs1%values = obs2%values
@@ -1541,7 +1626,17 @@
! Free up allocated storage in an observation type
type(obs_type), intent(inout) :: obs
-deallocate(obs%values, obs%qc)
+!write(*,*) 'destroy_obs: freeing space'
+if (associated(obs%values)) then
+ deallocate(obs%values)
+ nullify(obs%values)
+endif
+if (associated(obs%qc)) then
+ deallocate(obs%qc)
+ nullify(obs%qc)
+endif
+!if pointers are nullified() then this is safe (and simpler).
+!deallocate(obs%values, obs%qc)
call destroy_obs_def(obs%def)
end subroutine destroy_obs
More information about the Dart-dev
mailing list