[Dart-dev] [3797] DART/trunk/obs_sequence: In obs_sequence_mod, 2 important fixes.
nancy at ucar.edu
nancy at ucar.edu
Tue Mar 31 11:39:39 MDT 2009
An HTML attachment was scrubbed...
URL: http://mailman.ucar.edu/pipermail/dart-dev/attachments/20090331/e15244d4/attachment.html
-------------- next part --------------
Modified: DART/trunk/obs_sequence/obs_sequence_mod.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-03-31 15:32:06 UTC (rev 3796)
+++ DART/trunk/obs_sequence/obs_sequence_mod.f90 2009-03-31 17:39:38 UTC (rev 3797)
@@ -255,8 +255,10 @@
! Interactive creation of an observation sequence
type(obs_sequence_type) :: interactive_obs_sequence
-type(obs_type) :: obs
-integer :: max_num_obs, num_copies, num_qc, i, end_it_all
+type(obs_type) :: obs, prev_obs
+type(obs_def_type) :: obs_def
+type(time_type) :: obs_time, prev_time
+integer :: max_num_obs, num_copies, num_qc, i, end_it_all
write(*, *) 'Input upper bound on number of observations in sequence'
read(*, *) max_num_obs
@@ -282,8 +284,9 @@
! Initialize the obs variable
call init_obs(obs, num_copies, num_qc)
+call init_obs(prev_obs, num_copies, num_qc)
-! Loop to initialize each observation in turn; terminate by ???
+! Loop to initialize each observation in turn; terminate by -1
do i = 1, max_num_obs
write(*, *) 'input a -1 if there are no more obs'
read(*, *) end_it_all
@@ -293,11 +296,27 @@
if(i == 1) then
call insert_obs_in_seq(interactive_obs_sequence, obs)
else
- call insert_obs_in_seq(interactive_obs_sequence, obs, interactive_obs_sequence%obs(i - 1))
+ ! if this is not the first obs, make sure the time is larger
+ ! than the previous observation. if so, we can start the
+ ! linked list search at the location of the previous obs.
+ ! otherwise, we have to start at the beginning of the entire
+ ! sequence to be sure the obs are ordered correctly in
+ ! monotonically increasing times.
+ call get_obs_def(obs, obs_def)
+ obs_time = get_obs_def_time(obs_def)
+ call get_obs_def(prev_obs, obs_def)
+ prev_time = get_obs_def_time(obs_def)
+ if(prev_time > obs_time) then
+ call insert_obs_in_seq(interactive_obs_sequence, obs)
+ else
+ call insert_obs_in_seq(interactive_obs_sequence, obs, prev_obs)
+ endif
endif
+ prev_obs = obs
end do
call destroy_obs(obs)
+call destroy_obs(prev_obs)
end function interactive_obs_sequence
@@ -696,10 +715,35 @@
! Get the time for the observation
obs_time = get_obs_def_time(obs%def)
+! Assume we're starting at the beginning.
+! If we make this smarter eventually, here is where
+! we'd set the initial key number for a search.
+
+! If given an existing obs, be sure the new obs time is
+! consistent - later or equal to the given previous obs.
if(present(prev_obs)) then
- prev = prev_obs%key
- current = prev
- next = prev_obs%next_time
+ prev = prev_obs%key
+ current = prev
+ next = prev_obs%next_time
+
+ ! it is an error to try to insert an observation after an
+ ! existing obs which has a smaller timestamp.
+ if (prev /= -1) then
+ current_time = get_obs_def_time(seq%obs(prev)%def)
+ if (obs_time < current_time) then
+ !! or, do the insert searching from the start
+ !prev = -1
+ !current = -1
+ !next = seq%first_time
+ ! error out
+ write(msg_string,*) 'time of prev_obs cannot be > time of new obs'
+ call error_handler(E_ERR,'insert_obs_in_seq',msg_string,source,revision,revdate)
+ endif
+ endif
+
+ ! the insert code will search forward starting at the
+ ! given obs, so it is not an error to give an obs which
+ ! has a larger time than the next obs.
else
! Start search at beginning
prev = -1
@@ -784,7 +828,7 @@
obs_time = get_obs_def_time(obs%def)
last_time = get_obs_def_time(last_obs%def)
if(obs_time < last_time) then
- write(msg_string, *) 'tried to append an obs to sequence with bad time'
+ write(msg_string, *) 'time of appended obs cannot be < time of last obs in sequence'
call error_handler(E_ERR,'append_obs_to_seq',msg_string,source,revision,revdate)
endif
Modified: DART/trunk/obs_sequence/obs_sequence_tool.f90
===================================================================
--- DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-03-31 15:32:06 UTC (rev 3796)
+++ DART/trunk/obs_sequence/obs_sequence_tool.f90 2009-03-31 17:39:38 UTC (rev 3797)
@@ -331,6 +331,8 @@
call destroy_obs_sequence(seq_in)
if (remaining_obs_count == 0) then
process_file(i) = .false.
+ write(msgstring,*) 'No obs used from input sequence file ', trim(filename_seq(i))
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
cycle
else
process_file(i) = .true.
@@ -357,6 +359,9 @@
! pass 2:
+! blank line, start of actually creating output file
+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)
@@ -372,7 +377,7 @@
if (.not. process_file(i)) cycle
write(msgstring,*) 'Starting to process input sequence file ', trim(filename_seq(i))
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
call read_obs_seq(filename_seq(i), 0, 0, 0, seq_in)
@@ -412,10 +417,11 @@
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
- if (print_only) then
- call print_obs_seq(seq_in, filename_seq(i))
- endif
+ ! ensure the linked list times are in increasing time order
+ call validate_obs_seq_time(seq_in, filename_seq(i))
+ if (print_only) call print_obs_seq(seq_in, filename_seq(i))
+
!-------------------------------------------------------------
! Start to insert obs from sequence_in into sequence_out
!
@@ -467,8 +473,8 @@
total_num_inserted = total_num_inserted + num_inserted
else
- write(msgstring,*)'no first observation in ',trim(adjustl(filename_seq(i)))
- call error_handler(E_MSG,'obs_sequence_tool', msgstring, source, revision, revdate)
+ write(msgstring,*)'no first observation in ',trim(filename_seq(i))
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
endif
if (.not. print_only) then
@@ -484,7 +490,7 @@
enddo
write(msgstring,*) 'Starting to process output sequence file ', trim(filename_out)
-call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+call error_handler(E_MSG,'obs_sequence_tool',msgstring)
if (.not. print_only) then
print*, 'Total number of obs inserted : ', total_num_inserted
@@ -498,8 +504,7 @@
call write_obs_seq(seq_out, filename_out)
else
write(msgstring,*) 'Output sequence file not created; print_only in namelist is .true.'
- print *, trim(msgstring)
- call error_handler(E_MSG,'obs_sequence_tool', msgstring, source, revision, revdate)
+ call error_handler(E_MSG,'', msgstring)
endif
! Time to clean up
@@ -567,13 +572,13 @@
if ( num_copies1 /= num_copies2 ) then
write(msgstring2,*)'Different numbers of data copies found: ', &
num_copies1, ' vs ', num_copies2
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2, source, revision, revdate)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
num_copies = -1
endif
if ( num_qc1 /= num_qc2 ) then
write(msgstring2,*)'Different different numbers of QCs found: ', &
num_qc1, ' vs ', num_qc2
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2, source, revision, revdate)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
num_qc = -1
endif
if ( num_copies < 0 .or. num_qc < 0 ) then
@@ -586,12 +591,12 @@
if( str1 == str2 ) then
write(msgstring2,*)'metadata ',trim(adjustl(str1)), ' in both.'
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2, source, revision, revdate)
+ 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, source, revision, revdate)
+ 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, source, revision, revdate)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
endif
enddo MetaDataLoop
@@ -602,12 +607,12 @@
if( str1 == str2 ) then
write(msgstring2,*)'qc metadata ', trim(adjustl(str1)), ' in both.'
- call error_handler(E_MSG, 'obs_sequence_tool', msgstring2, source, revision, revdate)
+ 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, source, revision, revdate)
+ 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, source, revision, revdate)
+ call error_handler(E_MSG, 'obs_sequence_tool', msgstring2)
call error_handler(E_ERR, 'obs_sequence_tool', msgstring1, source, revision, revdate)
endif
enddo QCMetaData
@@ -637,7 +642,7 @@
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
' are before first_obs_days:first_obs_seconds'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -651,7 +656,7 @@
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
' are after last_obs_days:last_obs_seconds'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -671,7 +676,7 @@
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
' are on the discard obs_types list'
endif
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -688,13 +693,13 @@
else
!%! call select_obs_by_location(min_loc, max_loc, seq, all_gone)
msgstring = 'Select by general bounding box not implemented yet'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: no obs in ' // trim(seqfilename) // &
' are above the given height'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -716,13 +721,13 @@
if (.not. found) then
msgstring = 'QC metadata string: ' // trim(qc_metadata) // &
' not found in obs_seq file'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
' are outside the qc min/max range'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -743,13 +748,13 @@
if (.not. found) then
msgstring = 'Copy metadata string: ' // trim(copy_metadata) // &
' not found in obs_seq file'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
if(all_gone) then
if (print_msg) then
msgstring = 'Skipping: all obs in ' // trim(seqfilename) // &
' are outside the copy min/max range'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -763,7 +768,7 @@
if (print_msg) then
msgstring = 'Skipping: no obs in ' // trim(seqfilename) // &
' are above the GPS height threshold'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
endif
remaining_obs_count = 0
return
@@ -809,8 +814,8 @@
size_seq_in = get_num_obs(seq_in)
if (size_seq_in == 0) then
- msgstring = 'Obs_seq file is empty.'
- call error_handler(E_MSG,'obs_sequence_tool',msgstring,source,revision,revdate)
+ msgstring = 'Obs_seq file '//trim(filename)//' is empty.'
+ call error_handler(E_MSG,'obs_sequence_tool',msgstring)
return
endif
@@ -832,8 +837,8 @@
is_there_one = get_first_obs(seq_in, obs)
if ( .not. is_there_one ) then
- write(msgstring,*)'no first observation in ',trim(adjustl(filename))
- call error_handler(E_MSG,'obs_sequence_tool', msgstring, source, revision, revdate)
+ write(msgstring,*)'no first observation in ',trim(filename)
+ call error_handler(E_MSG,'obs_sequence_tool', msgstring)
endif
! process it here
@@ -898,6 +903,88 @@
end subroutine print_obs_seq
!---------------------------------------------------------------------
+subroutine validate_obs_seq_time(seq, filename)
+
+! this eventually belongs in the obs_seq_mod code, but for now
+! try it out here. we just fixed a hole in the interactive create
+! routine which would silently let you create out-of-time-order
+! linked lists, which gave no errors but didn't assimilate the
+! right obs at the right time when running filter. this runs
+! through the times in the entire sequence, ensuring they are
+! monotonically increasing in time. this should help catch any
+! bad files which were created with older versions of code.
+
+type(obs_sequence_type), intent(in) :: seq
+character(len=*), intent(in) :: filename
+
+type(obs_type) :: obs, next_obs
+type(obs_def_type) :: this_obs_def
+logical :: is_there_one, is_this_last
+integer :: size_seq, num_copies, num_qc
+integer :: i, j, key
+type(time_type) :: last_time, this_time
+
+
+! make sure there are obs left to process before going on.
+size_seq = get_num_obs(seq)
+if (size_seq == 0) then
+ msgstring = 'Obs_seq file '//trim(filename)//' is empty.'
+ call error_handler(E_MSG,'obs_sequence_tool:validate',msgstring)
+ return
+endif
+
+! Initialize individual observation variables
+call init_obs( obs, get_num_copies(seq), get_num_qc(seq))
+call init_obs(next_obs, get_num_copies(seq), get_num_qc(seq))
+
+!-------------------------------------------------------------
+! Start to process obs from seq
+!--------------------------------------------------------------
+is_there_one = get_first_obs(seq, obs)
+
+if ( .not. is_there_one ) then
+ write(msgstring,*)'no first observation in sequence ' // trim(filename)
+ call error_handler(E_MSG,'obs_sequence_tool:validate', msgstring, source, revision, revdate)
+endif
+
+call get_obs_def(obs, this_obs_def)
+last_time = get_obs_def_time(this_obs_def)
+
+
+is_this_last = .false.
+ObsLoop : do while ( .not. is_this_last)
+
+ call get_obs_def(obs, this_obs_def)
+ this_time = get_obs_def_time(this_obs_def)
+
+ if (last_time > this_time) then
+ ! bad time order of observations in linked list
+ call print_time(last_time, ' previous timestamp: ')
+ if (gregorian_cal) call print_date(last_time, ' Gregorian day: ')
+ call print_time(this_time, ' next timestamp: ')
+ if (gregorian_cal) call print_date(this_time, ' Gregorian day: ')
+
+ key = get_obs_key(obs)
+ write(msgstring,*)'obs number ', key, ' has earlier time than previous obs'
+ call error_handler(E_MSG,'obs_sequence_tool:validate', msgstring)
+ write(msgstring,*)'observations must be in increasing time order, file ' // trim(filename)
+ call error_handler(E_ERR,'obs_sequence_tool:validate', msgstring, source, revision, revdate)
+ endif
+
+ last_time = this_time
+
+ call get_next_obs(seq, obs, next_obs, is_this_last)
+ if (.not. is_this_last) obs = next_obs
+
+enddo ObsLoop
+
+! clean up
+call destroy_obs( obs)
+call destroy_obs(next_obs)
+
+end subroutine validate_obs_seq_time
+
+!---------------------------------------------------------------------
subroutine print_metadata(seq1, fname1)
!
More information about the Dart-dev
mailing list